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1. Introduction 

One often uses perturbation theory to produce predictions for the results of particle physics 
experiments in which the strong interaction is involved. Let us suppose that the measure¬ 
ment to be made is cast in the form of the cross section weighted by a function F(f ) that 
assigns a number to each possible final state /. The object then is to predict o[F\. Pertur¬ 
bation theory can be relevant if a high momentum transfer Q is involved in the reaction. 
In that case, one arranges the calculation so that one is performing an expansion in powers 
of a s (Q), which is small when Q is large even though a s (l GeV) is not small. If there 
are hadrons in the initial state, then low scale effects related to the initial state must be 
factored into parton distribution functions. Low scale effects related to the final state can 
be avoided if F is “infrared safe,” as described in Sec. | 2 ] and if its definition does not involve 
any small parameters. 

It will prove helpful to have an example, taken from e + e _ —► hadrons. Let F(f) be 
(1 — tf) 4 where tf is the thrust of the particles in /. If we call this observable T 4 , we have 

4 H] = . (1.1) 

A related example is obtained by collecting the particles in / into three jets using a suitable 
jet algorithm and defining my to be the mass of the most massive of the three jets. Then 
we can take F(f) to be dT^/dM 1 2 = (1 — tf) 4 5(m 2 — M 2 ), so that 

a[dT±/dM 2 ] = da[T 4 \/dM 2 = £dt (1 - t) 4 . (1.2) 

For fixed M 2 /s, this is an infrared safe observable, but if M 2 /s <C 1, its perturbative 
expansion contains large logarithms, log (M 2 /s). Quantities like this need a treatment 
beyond a fixed order of perturbation theory. In practical applications one often needs not 
only a treatment of the large logarithms but also a good model for how partons turn into 
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hadrons. That is, to make predictions for we need just short distance physics, but 

to understand da\T^\/dM 2 for all M 2 , we need to model also long distance physics. 

The most widely used tools available for perturbative calculations can be classified 
into two groups. The first is parton shower Monte Carlo (MC) event generators, such as 
HERWIG Q and PYTHIA [j2j], The second is next-to-leading order (NLO) programs such 
as NLOJET+- 1 - Q and MCFM 0. The NLO programs have NLO accuracy for hard cross 
sections like crpi], but do a bad job of predicting quantities like dcr[T±\/ dM 2 that represent 
the inner structure of jets. The parton shower Monte Carlo programs do not have NLO 
accuracy for hard cross sections but do a good job of predicting quantities like da\T^\/dAd 2 . 
Additionally, they have the advantage that their final state particles are hadrons rather 
than partons. There are also some NLO-MC hybrids. One example is the program of 
Frixione, Nason, and Webber i |0, which so far has been applied to cases with massless 
incoming partons but not to cases with massless final state partons at the Born level of 
calculation. The other example is that of ii®, which concerns three-jet observables 
in electron-positron annihilation and thus addresses massless final state partons but not 
massless initial state partons. This paper concerns a method for constructing NLO-MC 
hybrid programs. 

The shower Monte Carlo event generators can be made more powerful by incorporating 
the kT -jet matching scheme of Catani, Kuhn, Krauss and Webber 111 ]. This scheme is 


defined by first considering the cross section a\F\ to be divided into partial cross sections 

OO 

a[F\ = J2 °m[F] ■ (1.3) 


771=2 


Here a m [F} consists of the contribution to (t[F\ from final states that are classified as 
consisting of m jets according to a slightly modified version of the k? jet algorithm (as 
described in Sec. || of this paper). This jet algorithm depends on the choice of a resolution 
scale for the jets, which we call d- in \. 

A parton shower Monte Carlo program based on the &T-jet matching scheme is based 
on the Born squared matrix element for a m multiplied by a reweighting factor that we 
call W, which contains Sudakov factors of the sort found in parton shower Monte Carlo 
programs. This produces a realistic structure for jets at a scale above d- in \. The base 
calculation for a m [F] produces m partons. For each of these partons, the full program 
then produces a complete shower with d in i as a maximum scale. This fills in the details of 
jet structure at scales below d; n i. Both the factor W and the Sudakov factors that control 
the shower development at scales below d; n i depend on d- m i, but the product of these two 
factors is approximately independent of d- m i. 

Suppose that we use a program of this sort to calculate cr[F] for an infrared safe N- 
jet observable as defined in Sec. |2|. We assume that the definition F is such that the 
perturbative expansion of a[F] does not involve large logarithms. Thus F measures large 
momentum scale features of events with N or more jets. It vanishes if there are fewer than 
N narrow jets. As long as the jet resolution parameter di n i is appropriately chosen, the 
main contribution to cr[F] comes from a m [F] with m = N. This contribution is of order 
af rn , where B m = m —2 for e + e _ —► hadrons. For instance, in our example, T 4 is a three jet 
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observable. A calculation using PYTHIA with y/s = Mz shows that if we pick d; n i ~ 0.03, 
then 0-3 [T 4 ] gives the dominant contribution to crfTj, with <72 [T 4 ] ~ <74 [T 4 ] « 0.1 03 [T 4 ]. We 
discuss this further in Sec. ||. 

The perturbative expansion of a m [F] has the form 

cJm[F\ = C mfi [F] af m (Q) + C m> 1 [F] a^ + 1 (Q) + C m , 2 [F} a*™ + 2 (Q) + • • • . (1.4) 


(This applies for any m, not just m = N, but some of the coefficients are zero if m < 
N .) A leading order parton shower Monte Carlo program based on the kx- jet matching 
scheme will get the Born term, C m fi[F] af m (Q), in this expansion exactly right. Higher 
order contributions will be present, but getting the higher order perturbative contributions 
correctly is beyond the order of approximation intended in such a program. 

If we wanted more perturbative accuracy than described above, we might use an NLO 
calculation. There are, in fact, NLO calculations available for a wide variety of important 
processes. Modern calculations are in the form of computer programs, each of which is 
designed to work for a certain class of observables F. Such programs are very important 
for producing accurate predictions for the class of measurement functions F for which they 
were designed. Unfortunately ordinary NLO programs have significant limitations. For 
instance, an NLO program can calculate < 7 ( 14 ], but if one were to use the same program to 
calculate da[T,i]/dM 2 for M 2 / s -C 1, the result would not be even qualitatively right. For 
y/s = Mz, the physical distribution peaks at a jet mass of a few GeV and tends to zero for 
M —> 0. In contrast, the perturbative program gives for dcr\T/\/dM 2 a result that increases 
without bound as M 2 —> 0 and additionally contains a delta function at zero mass, 5(M 2 ), 


with a coefficient that is negative and infinite (see [ 10 ( for a similar example worked out 
numerically). 

The purpose of this paper is to extend to next-to-leading order the parton shower 
Monte Carlo idea based on the kx- jet matching scheme. That is, we want to keep the 
feature that each outgoing hard parton generates a full parton shower and hadronization. 
We also want to utilize the decomposition of a[F] into /2 a rn[F] so that we can include 
calculations for different numbers of hard jets in the same program. However, we now want 
the perturbative expansion of the calculated a m [F ] to be correct with respect to the first 
two coefficients in Eq. 0), namely C m $[F] and C m ^[F\. 1 

The general idea of the algorithm that we present applies, we believe, to lepton-lepton 
collisions, lepton-hadron collisions, and hadron-hadron collisions. This is because the dipole 
subtraction scheme that we use applies to all of these cases. However, we have elected to 
work out only the case of lepton-lepton collisions here and to leave the other cases for 
future publications. 

Before we launch into a construction that must be rather involved if it is to be precise, 
it may be useful to begin with an informal preview. Consider as an example the calculation 
a [T 4 ]. The dominant contribution to crfTjj] comes from <73 [X 4 ]. This quantity has a pertur¬ 
bative expansion of the form (1.4) with B 3 = 1. We seek a calculation that includes parton 


1 Of course, we can have NLO accuracy only for those values of m for which NLO calculations exist. For 

e + e“ —> hadrons, that would currently be m = 2, 3, and 4. 
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showering and hadronization and gets the two coefficients C’3,0 [T4] and 63,1 [X4] right. At 
the Born level ( a$ 3 ) there are three partons in the final state. At order a,f 3+1 there can 
be four partons and we might think of the four partons as arising from three partons by 
parton splitting, say q —■► q + g, even though the actual calculation involves the exact af 3+1 
matrix element. If we now add parton showering to one of the Born graphs we are in big 
trouble. The showering includes splittings q —> q + g, which we had already included. 
To include parton shower splittings of the Born level partons without double counting, we 
must expand their effect in perturbation theory and subtract the order af 3+1 contributions 
from the a^ 3+1 graphs. 

We see that there are subtractions from the order a,f 3+1 graphs that are derived from 
the splitting functions that generate the initial step of the parton shower. That is good, 
because NLO calculations inevitably involve subtractions. In this paper, we want to have an 
algorithm that can be used by NLO practitioners in a reasonably straightforward manner. 
For this reason, we design the parton splitting so that the subtractions are (with some 


minor changes) those of the dipole subtraction scheme of Catani and Seymour [ 12 1 . This 
scheme is expressed in a Lorentz covariant style and nicely expresses the complete available 
phase space for parton splitting. It is quite widely used for NLO calculations (for example 
in the programs NLOJET++ Q and MCFM Q). 

The one parton splitting just described represents the hardest splitting in the eventual 
shower. After this comes secondary showering that can be performed according to an ex¬ 
isting shower Monte Carlo program suitably modified to account for the starting condition 
that the hardest splitting has already occurred. 

Now, the main contribution to a calculation of erfTj] comes from final states that have 
three jets. The first splitting and the secondary showering generates the inner structure 
of these jets, which would be reflected in da\T^\jdM 1 . Experience has shown that quite a 
variety of showering algorithms are capable of doing a good job of representing this inner 
structure. In this paper, we leave open the choice of methods for the secondary showers. 
As long as the secondary showering has certain simple properties, the coefficients C 30 P 4 ] 
and ( 7 , 3,1 [Tj] will not be affected. 

With this structure, we hope to largely decouple the part of the calculation that needs 
NLO computations from the bulk of the parton showering and hadronization. We have in 
mind that NLO practitioners could then write code that could be used in conjunction with 
a variety of Monte Carlo event generators, even with Monte Carlo event generators that 
are written after the NLO code. 

We seek to be quite flexible not only with respect to the Monte Carlo event generator 
but also with respect to several of the function choices internal to the algorithm. The 
reason for this flexibility is to allow for improvements in the methods we suggest. 

The complete algorithm calculates a[F ] in the form 


1 N 

°[ F ] ~ Jj^2 W n F (.fn) 


(1.5) 


71=1 


That is, it is in the form of a Monte Carlo event generator with weights w n . As is natural 
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for a description that includes quantum interference, the weights for individual simulated 
events can be positive or negative. 2 

The algorithm presented in this paper combines two separate ideas. First, we use the 
kx -jet matching scheme to match calculations involving different numbers of hard jets [11]. 
One might call this m-jet/(m + l)-jet matching. Second, we use the dipole subtraction 
scheme 0 to match between the hardest splitting in a parton shower and the next-to- 
leading order contributions to the m-jet matrix elements. One might call this PS/NLO 
matching. Most of the paper concerns PS/NLO matching. The reader might want to 
work on this part first. For this purpose, one can simply take dim to be infinitesimally less 
than 1. In that case, the partial cross sections cr m [F] all vanish except for < 72 [F] and the 
reweighting factors W that appear in the paper can be taken to be IF = 1. Of course 
such a calculation would be of only limited phenomenological usefulness since it would 
apply only to two-jet infrared-safe observables F. To get to a method that applies to 
three-jet observables, we need to do something about the two-jet region and then we need 
m-jet/(m + l)-jet matching. 

In the following sections, we review the dipole algorithm [12] for doing NLO calcula¬ 
tions in the case of electron-positron annihilation. For readers who are familiar with this 
algorithm, our aim is to set up the notation that we use, which differs in some instances 
from that used in Ref. [0[. For other readers, our aim is to provide a compact intro¬ 
duction to the algorithm. We provide no proofs that the algorithm works. Furthermore, 
there are several functions that must be defined in order to fully specify the algorithm, 
but in the introductory sections that follow, we skip the formulas for these functions and 
merely state the important properties of the functions that follow from these formulas. 
The formulas are then summarized later on, in Secs. 00 and 0. Once the dipole sub¬ 
traction formalism has been set up, we define in Sec. |5] how to break up the cross section 
into partial cross sections involving different numbers of resolved jets, along the lines of 
Ref. 0- Then in Sec. |6| we modify the subtraction scheme just a little to accommodate 
this division. The main subject of the paper begins in Sec. [?] with an outline of the general 
structure for adding showers to the partial cross sections. This is followed by an exposition 
in Sec. ||, with several subsections, of how the showers are to be added while keeping track 
of the next-to-leading order correction terms needed to keep from changing the coefficient 
C m ,i[F\ in Eq. Q. The following sections contain details on the functions used in the 
dipole subtraction and showering constructions. We present some conclusions in Sec. 16. 


2. Notations 

In this section, we introduce some of the notations that we will use throughout the paper. 

We wish to describe the process e + + e~ —> hadrons. In a calculation at a finite order 
of perturbation theory, we consider final states consisting of n partons, with n > 2. We 
denote the momenta of these partons by {p} n = {pi,P 2 , ■ ■ ■ ,Pn}- We represent the phase 

2 In NLO programs, there are also weights w„, but the weights can be very large, with an event with a 
large positive weight being followed by a counter-event with a large negative weight. Here, the weights are 
not large. 
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space integration for the final state partons as 


(27t) 4 i5 4 

where Po is the momentum of the initial e + e~ pair, Po = (\/s, 0 , 0 , 0 ) if we use the c.m. 
frame. 

We denote the flavors of the partons by labels {f} n with /,; € {g, u, u, d, d,... }. Then 
the complete description of the final state momenta and flavors is specified by the list 
{p, f}n = We can define “addition” on flavors by saying 

that fi + f -2 = fz if there is a QCD vertex for /i + f -2 —*► fs- Thus, for instance, d + g = d 
and u + u = g, while d + u is not defined. The splitting of a final state parton with flavor 
fi can be represented by giving the pair of daughter flavors {fip, fi^} i n the set of flavor 
pairs such that f u + f U2 = fi • 

The matrix element A4 for a + b —► n partons depends on the spin and color indices 
of the n partons. In order to avoid writing these indices, we follow the notation of Catani 
and Seymour and write the matrix element as a vector 


^(Wn) = n((*0 3dj4pi S +frl)) 


1=1 


Pi - -Po 


( 2 . 1 ) 


| M({p,f} n )) 


( 2 . 2 ) 


in color < 8 > spin space. The squared matrix element, summed over colors and spins, is then 

(M({p, f}n)\M({p, f} n )) . (2.3) 

We consider a perturbative calculation that is designed to predict an infrared safe 
IV-jet observable. We need to be careful about what we mean by an infrared safe observ¬ 
able. There should be a function F n ({pi,p 2 , ■ ■ ■ ,p n }) defined for any number n of massless 
partons with momenta pi. The functions F n do not depend on the flavors of the partons 
and should be invariant under permutations of the momentum arguments pi,P2, ■ ■ ■ ,Pn- 
In a calculation at “all” orders of perturbation theory, the cross section for the observable 
would be 

a i F i=E s /«" ■ «• ifl d ° iPn ■ < 2 4 ) 

Here p® = \p t \ in F and da contains the delta functions for momentum and energy con¬ 
servation. The cross section needs to contain a regulator to control infrared divergences, 
which cancel between terms with different numbers n of partons. To ensure that these 
cancellations work, the F n functions for different values of n need to be related. Specif¬ 
ically, if two partons become collinear or one becomes soft (p t —> 0 ), the outcome of the 
measurement should be unaffected: 


Fn({Pl,P 2 ,-- • , (1 - \)Pn-lApn-l}) = F n _ 1 ({p 1 ,p 2 ,.. . ,p n - 1 }) (2.5) 

for 0 < A < 1. (The case A = 0 covers the possibility that parton n is soft and not 
necessarily collinear with any other parton.) Cancellation of infrared divergences also 
requires that the measurement is unaffected when many partons become soft or group 
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themselves into collinear groups (jets). To cover this possibility, we also require that the 
functions F n ({p} n ) be smooth functions of their arguments. 3 

Ultimately, one wants to apply the measurement function to hadrons. Indeed, we want 
to do that in this paper after supplementing the perturbative calculation with a simulation 
of parton showering and hadronization. For this purpose, one uses the same functions F n , 
but now applied to hadron momenta pj with pj > 0. The extensions of the F n to pj > 0 
should then have smooth limits as any or all of the pj approach 0. 

We have defined what we mean by an infrared safe observable. In this paper, we 
consider an infrared-safe “IV-jet” observable. This means that additionally 


Fm({Pi,P2, ■ ■ ■ ,Pm}) = 0 , m < N . (2.6) 

Throughout the paper we use the standard notations Cf = (Nj — 1)/2N C , C\ = N c , 
Tp; = 1/2 for color SU(IV C ). We assume that there are Nf flavors ( + Nf antiflavors) of 
quarks in the fundamental representation. We typically use the notation / to represent a 
parton flavor, / G {g, u, u, d,... }. Then we use coefficients Cf, Kf, and 7 / defined by 


Cf 

C g 

7/ 


7g 

K f 


Cf for / G {u, u, d,... } , 

C A , 

§ C F for / G {u, u, d,...} , 

TtCa - -T R N { , 

0 6 

7 7 r 2 \ 

- - y j C F for / G {u, u, d,... } 


67 7T 2 


10 , 


, C A - T r N{ . 

18 6 / 9 


(2.7) 


The running coupling is a s (p) evaluated at scale p, often a transverse momentum. 
When no scale is indicated, we mean a s = a s (p R ) where the p R is a Hxed renormalization 
scale, usually chosen as some fraction of y/s. For our NLO calculations, we would use the 
two loop running coupling. In construction our matching with parton showers, we will 
require the first order relation 


where /3q = 7 g . 


«s(h) 

a s (m) 


1 - Po log 



« s (m) 

2 vr 


+ ^(“s) 


( 2 . 8 ) 


3. Construction and deconstruction of parton splitting 


The dipole algorithm of Catani and Seymour [12] is based on a physical picture involving 
parton splitting, which turns m partons into m + 1 partons. Deconstructing the splitting 


3 A smooth function has an infinite number of derivatives, and thus a Taylor expansion to any order. 
This property is stronger than one really needs, but successively weaker requirements become successively 
more unwieldy. 






turns the m + 1 partons back into m partons. In this section, we describe the splitting 
kinematics without giving the detailed formulas from Ref. [ 12], which are given in Sec. 10. 

We describe deconstruction first. Suppose that we have a list of m +1 parton momenta 
and flavors, {p, f} m +i • One imagines that partons i and j are produced by the splitting of 
a mother parton with flavor fij and momentum pij. We need one more parton, with index 
k, to describe the splitting in the scheme of Catani and Seymour. Parton A; is a “spectator 
parton” that absorbs some momentum associated with the splitting. 4 

Consider first the flavors. The mother parton has flavor fij = fi + fj . The flavor 
of the spectator parton is not changed: /& = //-. There is less information in the list of 
just one flavor fij than in the list of two flavors The missing information is the 

flavor splitting choice, which can be specified by giving the pair in the set with 

fi + fj = fij ■ 

We now extend this idea to the momenta. The massless momenta {pi,Pj,Pk} determine 
new momenta {pij,Pk} of just two on-shell massless partons together with three splitting 
variables. The structure of this transformation is simple in the limiting case of collinear 
splitting. If Pi- f>j = 0, we have 


Pij =Pi+ Pj (when pi ■ pj = 0) , 
while the momentum of the spectator parton remains unchanged, 


(3.1) 


Pk = Pk (when fi ■ pj = 0) 


(3.2) 


It is not possible to retain these relations away from the collinear limit. However, Catani 
and Seymour still maintain 

Pij +Pk=Pi+ Pj + Pk (3.3) 

while keeping all of the momenta massless. 

Since one eliminates three degrees of freedom in going from {pi,Pj,Pk} to {pij,Pk}, we 
can supplement {pij,Pk} with three splitting variables. For our purposes, it is convenient to 
call these y, z , <p. The most important role in the formalism is played by the dimensionless 
variable y proportional to the virtuality of the splitting, so that 


y = 0 if and only if pi ■ pj = 0 . (3.4) 

The variable z is a momentum fraction representing the share of the mother parton mo¬ 
mentum that is carried by parton i and <f is an azimuthal angle. 

The discussion above can be summarized by saying that there is a map 

{(.PiJi), ( PjJj), {PkJk)} -> {(PijJij), ( PkJk),y,z,fJiJj} • (3.5) 

Parton splitting is the other way around. Here we begin with a list of the momenta 
and flavors {p, f} m of m partons. We imagine that one of these, parton l, splits, producing 

4 The spectator parton also plays a role in the color flow associated with the splitting, as expressed in 
the color matrices in the definition of T>ij,k below. 
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daughter partons with with flavors and momenta ■Pvi}- We again need 

a spectator parton, with index k. To specify the splitting we need splitting parameters 
{y, z, 4>, fi t i, fi t 2 }, with // i 1 + fi 2 = fi ■ Then Catani and Seymour specify a map 


{(Pi, fl), ( Pk , /fe), V, z, 4>, fi t 1 , /z, 2 } -»• /z,i), (pz, 2 , fl, 2 ), (Pfc, /fe)} 


(3.6) 


The map ( |3.6| ) is precisely the inverse of the map (3.5). Only the variable names are 
different. The part of these maps that concerns the flavor splitting is trivial. The part that 


concerns the momentum splitting is not trivial. The formulas from Ref. [12] are given in 


Sec. 10 


These maps can be rather trivially extended to include all of the partons. Suppose 
that we start with a list {p, f} m of m parton momenta and flavors and that we want to 
split parton l with the help of spectator parton k using splitting variables {y, z, <j >,,/),] , fi^}- 
We need /).] + /;. 2 = fi ■ We can produce a list {p, f} m +i of m + 1 parton momenta and 
flavors by removing partons l and k from the original list and adding p/ ; 1 , /; > 1 , p^ 2 , fl, 2 and 
pk, fk from Eq. ( |3.6| ) to the end of the list. Then {p, f} m+ 1 equals the new list. It will 
prove useful to call the complete transformation 7Z^k, 


{.P, f }m+i = T^i,k{{p, f}m, {y, z , <t>, fa 1 , fi )2 }) ■ 


(3.7) 


Now suppose that we start with a list {p, f} m + 1 of m + 1 parton momenta and flavors 
and that we want to combine partons i and j with the help of spectator parton k. We can 
produce a list {p, f} m of m parton momenta and flavors by removing partons i,j and k 


from the original list and adding Pij,fij and pk,fk from Eq. (3J5) to the end of the list. 
Then {p, f} m equals the new list. We also get the splitting variables {y, z, (j), fi, fj} with 
the help of Eq. (|3.5|). It will prove useful to call the complete transformation Qij.k, 


{{P, f }m, {y, Z, 4>, fi, fj}} = Qij,k({P, f}m+ 1 ) 


(3.8) 


The transformations 1Z and Q are inverses of each other in the sense that if we supply 
the right permutation operators II (i,j, k ) and II(Z, k') we have II (i,j, k)TZi : k'^{l, k')Q r] ,k = 
1. The permutations will not be of much concern to us since the functions that we use that 
are functions of {p, /} are invariant under permutations of the parton labels. 

Where possible, for a final state parton l, we denote the complete set of splitting 
variables by a single letter, 

Y i = {yi,zi,<i>i,fi,i, fi,2} ■ (3.9) 

With this notation we can abbreviate 



d(fi 1 j/; 

27t 2 z—fi,i+fi,2 



(3.10) 


Now that we understand the parton splitting, we are ready to examine the construction 
of the cross section at NLO. 
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4. Structure of the NLO cross section 


We consider an iV-jet cross section in electron-positron collisions. The cross section correct 
to NLO is based on a tree level cross section for the production of m partons, with m = N, 
together with one loop graphs for the production of m partons and tree graphs for the 
production of m + 1 partons. The cross section is constructed as a sum of terms 

(J n l° = (J b + c7 r- a + (T v + a ( 4 .1) 

In the first term there is an integration over m-parton phase space, while in the second 
term there is an integration over m + 1-parton phase space. The first term is the Born 
contribution, proportional to = a™ -2 . The second term is a correction proportional 
to af m+l associated with real parton emission, which comes with a subtraction term that 
eliminates the soft and collinear divergences. The third term is a correction proportional 
to af m+1 that is associated with a virtual parton loop. There are m partons in the final 
state. There is a corresponding subtraction term in which there is an integration over 
the phase space for one parton, which is performed analytically in 4 — 2e dimensions to 
produce 1/e 2 and 1/e terms that cancel 1/e 2 and 1/e terms that would be present without 
the subtraction. 

We will begin with the Born contribution. 


4.1 The Born contribution 

The Born contribution takes the form 

= ^ / 'dr({p} m )(M({p,f}m)\M({p,f}m))F m ({p } m ) . (4.2) 

^ {f}m J 

There is an integration over the final state momenta and a sum over final state flavors {f} m 
with a symmetry factor 1/m!. Next is the squared matrix element for the production of 
the m final state partons. Finally, there is a final state measurement function, F m ({p} m ). 


4.2 The real emission contribution 

The real emission contribution J dcr R along with its subtraction f da A has the form 


<7 


R-A 


/ m+1 


[da R — da A 


(mTTjT 51 /<«WW.) 

{f}m +1 

f}m+l)\M({p, f}m+1 )> Fm+1 ({]5} m+1 ) 


-EE 

i,j k^i,j 

pairs 


(4.3) 


The real emission term der R is represented by the first term in square brackets. Here we have 
the same sums and integrals as in cr B except that now there are m + 1 final state partons. 
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The momentum and flavor variables are all written with hats in order to distinguish them 
from the m-parton variables. There is the squared amplitude to produce the m + 1 partons. 
The resulting partonic cross section is multiplied by the measurement function F for m + 1 
partons. 

The contribution da R has a potential singularity when any of the dot products pi ■ pj 
tends to zero. In the second term in the square brackets, representing dcr A , we sum over 
subtractions designed to reduce the strength of these singularities. The subtractions are 
labeled by (unordered) pairs {i. j} where i.j £ {1 ,m + 1}. For each choice of {i,j}, 
there are a number of subtraction terms labeled by the index A: of a spectator parton, which 
can be any of the partons except i or j. We use the splitting deconstruction transformation 
from Eq. (|3.8| ) to define m parton variables and splitting variables according to 

{{ P , m k , Y ijik } = Qt jt k({p, f}rn+1 ) . (4.4) 


The ij, k superscripts or subscripts remind us that the definition of these variables is 
different for each choice of parton indices i,j,k. 

The subtraction term in Eq. (|4.3| ) contains a dipole subtraction function We 

state the definitions in Sec. 11, Eqs. (11 .2 . . B |11.4 |11.5| , |11.6| , |l 1 . 7] ) , but here simply 
note that has the structure 


Pij,fc({P, /}m+l) — 


2 Pi ■ Pj 




Tjj • Tfc 

_ r p2 

ij 


Vij{{pij,y, z, 4>}ij,k, fi, fj) \M({p, f}% k )) 


(4.5) 


The function is based on the Born level amplitude |A 4({p, f}m k )) for the flavors and 
momenta with partons i and j combined. There is an operator T t j ■ T k /T ri 2 j that acts on 
the color degrees of freedom of \M). Specifically, T"- is the infinitesimal SU(3) generator 
matrix in the a direction acting on the color indices for the mother parton that results 
from combining partons i and j (which, to be precise, has been placed in the second- 
to-last parton slot by the transformation Qij^k)- Similarly, TP is the generator matrix 
acting on the color indices for the spectator parton. Then ■ T k = ^) a T)“T^. In the 
denominator, Tf- = is a number, Cf or C\, depending on whether parton ij is 

quark or antiquark or a gluon. There is an important identity involving the color matrices. 
Invariance of \M) under color rotations implies that (Ty + T k )\PA) = 0. Thus 


Tij • T k 

2^—tT 

k^ij W 


= 1 


(4.6) 


when operating on |A4). Next in Eq. (|4.5P , there is an operator Mj.k that depends on 
the momenta and acts on the spin degrees of freedom of \M). The dipole subtraction 
functions 'Dij k are constructed so that their sum matches the squared matrix element 
{M{{p, i n the limit in which the matrix element is singular. 

Taking some liberties with the notation, we will also use P;,fc({p, f}m] Y) to denote 
the same function written in terms of the variables {{p, /} m ; T} given by Eq. (|3.8|). 
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Next in da A there is a measurement function Fj{{p}m k ) evaluated at the momenta for 
the m-parton state. The two measurement functions in the two terms in Eq. ( |4.3| ) match 
in the limits in which the matrix element is singular because of the infrared safety property 

©■ 

4.3 The virtual loop contribution 

The virtual loop contribution along with its counterterm has the form 

cr V+A = / [da v + f da A ] = ^ / dr{{p} m ) V({p,f} m ) F m ({p} m ) . (4.7) 

J m J 1 J 

XJ sm 

The function V comes from the one loop matrix element. The matrix element has been 
calculated analytically in 4 — 2e dimensions. Then its l/e n pole terms and certain finite 
pieces have been subtracted. The result V is constructed from Born amplitude and certain 
spin and flavor dependent functions. The terms subtracted are precisely what was sub¬ 
tracted from the real emission contribution, but with the opposite sign. The details of V 
are laid out in Sec. [T^, Eqs. (12 .2 , 12.3 , 12.4 ). 

5. Partial cross sections 

In subsequent sections, we will add parton showers to our calculation. Before we do this, 
however, it is useful to divide the cross section into partial cross sections a m with m = 
2,3,... that are based on the cross sections for 2 partons —► m partons and have higher 
order corrections added. This construction is based on the kx -jet matching scheme of 
Catani, Krauss, Kuhn, and Webber 01- 

We consider the calculation of an infrared safe iV-jet observable, as described in Sec. ||]. 
In particular, the measurement function F obeys Eq. (^Ij), F m ({pi,p 2 , ■ ■ ■ ,p m }) = 0 for 
m < N. We organize this calculation according to the number of partons in the final state, 



Here G n has a lowest order contribution proportional to a ^ n , where B n = n — 2, 

G° n = (M({p,f} n )\M({p,f} n )) . (5.2) 

There are higher order contributions, proportional to af n+L with L virtual loops. When 
the cross section is expanded in this form, some kind of regulation is needed on the integrals; 
the divergences then cancel between terms with different numbers of partons n. 

Given an n-parton final state, {p, f} n , we group the partons into jets. We use the 
recursive jet finding algorithm Jl3|, modified slightly to make use of the parton flavor 
information available. At each stage of the algorithm there is a dimensionless jet resolution 
function dij that approximates (for small angles and virtualities) the squared transverse 


(12.2, 12.3, 12.4 
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momentum of one jet (group of partons or a single parton) with label i with respect to 
another with label j with which it might be grouped, divided by s: 


dij = d(pi,pj) 


^ min [E?, E?] (1 - cos (%)) 


(5.3) 


Here we use the e + e 


c.m. frame to define 6 Q as the angle between p t and pj and to define 


Ei = p^, Ej = p°j. 

We begin with n partons. At each stage of the algorithm, we combine two partons, so 
that as the algorithm progresses, there are fewer and fewer partons left. There are some 
rules. In order to combine partons with flavors ft and / c , there must be a flavor f a such that 
there is a QCD vertex for f a —> fife- In addition, the initial splitting in e + e - —► hadrons 
must be 7 jZ —> qq. Therefore, when working backwards from the final state, it is not 
allowed to combine the last remaining q-q pair into a gluon. Thus at each stage, we find 
all pairs {i,j} of partons that are allowed to be combined by these flavor considerations. 
For each such pair, we calculate the distance measure d. tJ given in Eq. ( |5.3|) . The allowed 
pair with the smallest dij is combined by adding the four-momenta of the daughter partons 
to form the four-momentum of the mother and determining the flavor of the mother such 
that there is a QCD vertex for the mother to split to form the daughters. 

This process gives a sequence of resolution parameters dj({p, f} n ) at which two jets 
were joined, reducing J jets to J — 1 jets. Typically dj < d j_ 1 . However, the flavor consid¬ 
erations discussed above may invalidate this ordering when we must reject the parton pair 
with the smallest dij and choose a pair to combine that has a larger dij. This misorder- 
ing does not happen in the leading approximation in a parton shower. If dj_i calculated 
according to Eq. (|T^) is smaller than dj, we simply redefine dj- 1 to equal dj. Then the 
definition gives d n < d n _i < ■ ■ ■ < d 3 . We also define dj({p, f} n ) = 0 for J > n. 

We can use this kx jet algorithm to divide the cross section into partial cross sections 
<7 m with specified integration ranges, 


a[F\ = J2 °m[E] , (5-4) 

m=2 


with 

00 1 r 

CTmlF] = dr ({p}n) 0(d m+1 ({p, f} n ) < d] n ] < d m {{p,f}n )) 

n ' {fU J 

X G n ({p, /}„) F n ({p} n ) . 


(5.5) 


In G m there are exactly m jets that are resolvable at a scale dim that can be chosen to suit 
our purposes, as discussed briefly in Sec. |I]. Although there may be more than m partons, 
there are not more resolvable jets. There are no infrared divergences in a m arising from 
two of the m jets becoming collinear or one of them becoming soft because the singular 
region is removed by the cut dim < d m ({p, f} n )- There are also no infrared divergences in 
a m arising from the possible subjets becoming collinear or soft because of the cancellation 
between real and virtual graphs. 
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Assume that F is an infrared-safe iV-jet observable that is sensitive only to event 
structure at a large momentum scale Q ~ y/s. Then as long as the jet resolution parameter 
dini is appropriately chosen, the main contribution to a[F] comes from ctn[F], with N jets 
resolvable at scale d\ n \. We touched on this topic in Sec. || and discuss it further in Sec. [14]. 

Our aim in this paper will be to add showers to the perturbative calculation of a m [F] in 
such a way that the first two terms in the perturbative expansion of <J m [F] are reproduced 
without “double counting” between the shower splittings and the splittings that are part 
of the contributions to (j m \F] of order First, however, we need to indicate how the 

dipole subtraction scheme for NLO perturbative calculations can be modified to work with 
the partial cross sections a m [F], 


6 . Partial cross sections with dipole subtractions 


The integrals in Eq. ( |5.5| ) still need regulation because the real-virtual cancellations happen 
between terms with different values of n. However, as long as we limit the calculation to 
next-to-leading order, it is simple to adapt the dipole subtraction scheme to make the 
cancellations happen inside of each of the integrations. The idea is to write each partial 
cross section a m as a sum of three terms, with corrections suppressed by two powers of a s , 

°m = <4 + 0'm"' A + a m +A + 0(a^ m+2 ) . (6.1) 


For the Born contribution, we simply modify Eq. (|4.2|) by inserting the appropriate 


cut, 


dr({p} m )0(dini < dm{{p,f}m )) 

m\ z ' / 


{/}« 


( 6 . 2 ) 


X f}m)\M{{p, f} m )) F m ({p} m ) 


(Note that there is no cut d m+ \ < di n i since there are only m partons.) 

In the real emission contribution, we follow Eq. (|4.3|) , inserting the d m+ \ < di n i < d m 
cut in the main term and a similar cut in the subtraction term, 


<J 


R-A 


(m + 


£ jdrm mH ) 


{f}m+1 

(M({pJ} m +i)\M({p, f} m +i)) F m+ i({(p} m+ i) 

X ({)P, /}m+l) ^ ^ini ^ /}m+l)) 

-EE 

hi k^i,j 

pairs 

x 9(d({p, f}% k ,{l,y,z}ij,k) < d ini < d m {{pj }% k )) 


(6.3) 


In the subtraction term, there is a sum over parton pairs {ij} and spectator partons k. 
For each term, there is an appropriate cut on the momenta. First, the m-parton state 
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that results from combining partons i and j must be resolvable at scale d\ n j. Second, the 
splitting of the mother parton thus obtained must be unresolvable according to a resolution 
function d, 


d({p, f} m , l,yi,zi) = —yi min 
s 


1 - Zl 
Zl 


Zl 


1 - Zl 


(6.4) 


Here si is a virtuality scale appropriate to parton l. The simplest choice would be si = s. 
However, for our purposes, it will prove useful to choose a value obtained from {p, f} m 


that we will specify in Sec. 13, Eq. (13.16). Aside from the factor si/s, the function 
d({p, f} m J,yi, z i) is derived from the resolution variable of the krp jet finding algorithm, 
Eq.(^lJ). Except for a y and z independent factor, it gives an approximate version of 
the resolution variable associated with the splitting {p} m —> {p } m +1 generated with the 
splitting parameters yi , zi with the help of spectator parton ki , 


d(pi,uPi, 2 ) = 2pi Pkl d({p,f} m ,l,yi,zi) [i + <D{y/m)\ 


Si 


(6.5) 


We describe this more precisely in Sec. 10 at Eq. ( 10.19|) . The effect of the cuts on the 
dipole terms is easiest to understand in the case that for the m + 1-parton state that we 
start with, one pair has a very small resolution parameter, while once that pair is combined 
the other pairs are well separated. In this case, the cuts provide that only the dipole term 
for the pair with the small resolution parameter contributes. Note that the cancellation 
needed as any two of the partons {p} m - i-i become collinear or one of them becomes soft is 
left intact. 

We mention here a subtle point with respect to infrared singularities in Eq. (|6.3|). 
Consider the {ij. k\ subtraction term. As we integrate over {p. /} m +i, we can encounter 
a point at which two of the momenta in {p]m are collinear or one is zero. (The simplest 
situation is that two of the original {p} m+ i other than pi and pj are collinear or one is zero.) 
At this point, Aj.fc oc {M({p, f}m k )\M.{{p, f}m k )} is singular. With an IV-jet observable 
and with m = N , the function F({p} l Jj k ) vanishes at this point, so the singularity is 
weakened or eliminated. This is what happens in the standard dipole subtraction method 
im . However, suppose that m > N. Then F({p}fy k ) does not vanish at the point in 
question, but still d m ({p, f}m) vanishes, so that the cut di n i < d m ({p, f}m ) eliminates 
the singularity. 

In the virtual contribution, we write 


cr, 


V+A _ 


1 


ml 


/ dr({p } m ) 0 (d ini < d m ({p,f}m)) F m ({p} r 


{fh 


V({p, f}m ) - EE Cl, k ({ P J}m,d h 

l k^l 


( 6 . 6 ) 


This follows Eq. (0). with the addition of a cut d - m ; < d m . There are also added terms 
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involving functions Cj_ ^ 


Ci,k({pJ}m,d ini ) = f — f dz f < d({p,f} m ,l,y,z)) 

° V J ° ^ h,h 1 2 (6.7) 

x y( 1 -J d J 2 PrPk v lA{p ju Y ) . 

This arises because we modified the dipole subtractions in Eq. ( |6.3D by imposing an extra 
cut d({p, /}nf fc , {l, y, z}ij t k) < d[ n [ that was not present in Eq. (4.3). What we removed from 


the real emission subtraction we must also remove from the virtual loop subtraction. Doing 
that and changing integration variables from {p, f } m +1 to ({p, f} m ; Y{) gives Eq. ( |6.6|) . The 
details of the variable transformation can be found in Sec. 0, Eq. (|10.17|) . The functions 
Ci,k{{Pi f}m, dim) are evaluated in Sec. |l^, Eqs. ( 12.8 , 12. 9| ) 

Using dipole subtraction with cuts as defined above, one could construct a computer 
code that would calculate at next-to-leading order the expectation values of infrared safe 
IV-jet measurement functions for N = 2,3,4,... up to the value for which one had the 
required calculated matrix elements. The program would produce partonic events with 
weights, and the same events could be used for the calculation of different IV-jet observables 
with different values of N. The practical value of such a program would be minimal, since 
it would not add anything to having separate programs for each value of N. However, if 
one could add parton showers to the calculation thus organized, the program would have 
some added value over the separate perturbative programs. It is to this goal that we now 
turn. 


7. Partial cross sections with showers 


In the following sections, we discuss the construction of a parton shower algorithm that 
matches the dipole subtraction scheme for NLO calculations. We incorporate the A^-jet 


matching scheme [11], so that the cross section is generated as a sum of partial cross 
sections a m as in Secs. || and |6], but this time with showers added. We suppose that we 
can calculate the perturbative m-parton cross section, <r^, at the Born level, af™ = 
for m = 2,3,..., m max . Furthermore, we suppose that we can calculate the order af m+l 
corrections, cr^ l ~ A and <r^ +A , for m = 2,3,..., tunlo- The values of uinlo and m max 
depend on our ability to calculate, but we can always assume that m max > tunlo + 1 - 
The cross section computed with parton showers will consist of contributions from each 
available m, 


771NLO 777. m ax 

a NLO+S = g [ a B+S + a R+S + a V+S] + £ a B+S (7. 1} 

m =2 ra=mNLo+l 

For the contributions at NLO level, there are three terms, which correspond to Born, 
real emission, and virtual loop terms with showers added (“+ S”). For the remaining 
contributions there is only a Born term. We will arrange that (for a suitably behaved 
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observable) 


,B+S 

J m 


r R+S , _V+S 

' u m 


_ , _R-A , _V+A 

®m T T ® m 

+ 0(a*™ +2 ) + 0[af' 


x 1 GeV/ yfs ) 


(7.2) 


and 


(jB+s = a B + 0( a B m +1) + x ]_ GeV/Vi) 


(7.3) 


The main body of this construction is in the following section (with several subsec¬ 
tions). There we construct (7® +s , cr^ +s , and u^ +s for 2 < m < tunlo and show that these 
quantities sum to a^ + o-^~ A + cj^ +a to NLO accuracy. We also obtain the leading-order 


result in Eq. (7.2) as a byproduct. 


8. The shower construction 


In this section, we discuss <tB+ s , cjB+ s , and (7^ +s f° r 2 < m < tbnlo and show that these 
quantities sum to + <7 ^ _A + ct^ +a to NLO accuracy. We begin with <7® +s in Sec. S.l 
and the following subsections. Then <7^ +s is described in Sec. |8.7|. The virtual contribution 


cjV +a is described in Sec. 


3.8. We sum the contributions in Sec. 8.9 


8.1 Born term with showers 

Our discussion begins in this subsection with erB+ s . We define 



x I{{p,f}m',l,k,Yi) 


This formula is illustrated in Fig. |l]. The first line contains integrals over Born level parton 
momenta and a corresponding sum over parton flavors. The second line contains sums over 
choices l of the parton that splits and a spectator parton k along with an integral over the 
splitting variables 1) and a matrix element of certain operators Ei.k acting on the Born 
amplitude \M). The operators Ej ^ together with the other factors in the formula describe 
the formation of showers from the Born level partons. 

Our goal in the following subsections is first to explain the parts of Eq. Q and 
then to show that, with suitable definitions of the functions in Eq. o, the perturbative 
expansion of (jB+ s gives the Born cross section er® at order a ^ rn . Expanding to one more 
order, we get certain order af m+1 terms that we need to keep track of. We will then be 
able to add correction terms that allow us to match the contributions in the dipole 

subtraction scheme. 
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Figure 1: Illustration of cr® +s , Eq. (8T). The basic Born graph lies outside the dashed lines. 
The three partons are required to be resolvable at scale G?j n i. To the left of the first dashed line, 
we illustrate one contribution to | M). To the right, we illustrate one contribution to (_Ad|. We 
multiply by a reweighting factor W, not illustrated, that contains Sudakov factors for a constructed 
parton splitting history. Between the dashed lines there is a splitting of one of the partons l, using a 
spectator parton fc, as described by the splitting operator Ei^. Finally, there is secondary showering 


represented in the figure by diamonds and in Eq. (8T) by the function I. The diamonds split the 
partons into a shower of partons and then into jets of hadrons. 


8.2 The Born level partons 


The first line of Eq. (8.1) concerns what we call the Born level partons. These are labelled 
by an index l that takes the values l G {1, 2,..., m}. The Born level partons are produced 
in the hard scattering and evolve into final state showers. 

There is an integration over the phase space for the m final state partons with momenta 
{p}m- There is also a sum over the flavors {f} m of the final state partons. The integration 
over the momenta {p} m is restricted by a factor 

d(d; n i ^ d m ({p, /}m)) • 

Here d m ({p, f} m ) is defined by applying the kx jet finding algorithm to the m parton 
momenta, as in Eq. ( [~>.3p . As discussed in Secs. || and |], we classify a final state of n partons 
by the values of the jet resolution parameters d m . The contribution to the observable from 
states with d m +i < dim < d m are calculated using cr^ +s and its corrections. The cut on 
d m , in Eq. (8T) implements the first part of this general plan. 


There is also a factor W m ({p, f} m ), which is the product of factors associated with the 


splitting history that matches the found jet structure, following the method of Ref. [11]. 
Imagine that the found jet structure had been generated by parton shower Monte Carlo 
program. Then for each vertex V in the graph that represents the splitting history, there 
would be a factor a s (\/dvs), where dy is the kx scale of the splitting at that vertex. In 
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the exact squared matrix element (J\A({p, /} m )|.A4({p, /} m )), there is a factor a s (/hi) for 
each vertex. We switch from a s (//R.) to a s (y/dys) by including in W m ({p, f} m ) a factor 


a s (\Zdys) 

a s (m) 


( 8 . 2 ) 


for each vertex. The parton shower Monte Carlo program would also generate a Sudakov 
factor that represents the probability that the partons did not split between each scale dy 
at which a splitting occurred and the next smaller scale dy and the probability that the 
partons did not split between the smallest scale dy at which a splitting occurred and the 
limiting scale di n i. This Sudakov factor is also included in W m ({p, f} m ). The details are 
given in Sec. 13. 

All that we need to know at present about W m ({p, f} r 
expansion in powers of a s (pn) so that 


is that it has a perturbative 


W m ({pJ} m ) = l + ^^w£\{p,f} m ) + --- . (8.3) 

We will use this property in our perturbative analysis. 

In the second line of Eq. ( |S.lD there is a matrix element of certain operators E^k that 
act on the Born amplitude | M({p, f} m )) to generate the splitting of the partons. Our next 
task is to describe this parton splitting. 


8.3 Description of parton splitting 

Each Born level parton has the opportunity to split into two daughter partons. The one 
that splits with the largest (suitably defined) evolution parameter is designated as parton l. 
The splitting of the remaining partons is left to be described by the function / in Eq. (|3.1|). 
This splitting of parton l with the help of spectator parton k is described by an operator 
Ei t k- 

There are a number of parameters that describe the splitting of parton l, 


Yi = {yi,z h (j>i, 2 } • 


(8.4) 


as in Eq. 


We integrate over these variables in Eq. (0), using f dYi as defined 


m 


Eq. Mil) . 

The momenta and flavors of the m + 1 daughter partons are given in terms of the 
momenta and flavors of the m mother partons and the splitting variables Y) by the trans¬ 
formation ( 0 ), 

{p,f}m+i = Ki,k({p, f}m,Yi) . (8.5) 


8.4 The splitting functions 

Now, we turn to the splitting function £ 7 ^ in Eq. (|S.lD, which is an operator on the color 


- 20 








and spin space of parton l in the vector \M). This operator has the following form, 


Eik = 


/»oo 

/ drd(r - Ri{{p, f} m , yi, zi)) 0(d({p, f} m , l, y t , z t ) < d h 
Jo 

Ti ■ T k a s (k T ({p, f} m , l, V, z)) 


rp 2 

' ± l 


2t r 

r-l i / /• 1 




x exp| — J dr'J J^dz'^5{r' - R v {{p, f} m ,y',z')) 


( 8 . 6 ) 


v 


X /}m,i ,y ,z) <du 

X < S{y',z',fi /)) 


«s(Mfa/}rn/,?A^)) 

2 vr 


The parton splitting is organized according to an evolution parameter r, which is defined 
in Eq. ( |8.6| ) to be a certain function i?/ of the hard parton momenta and flavors { p , /} m 
and the splitting parameters yi, Z[. The function Ri measures the hardness of the splitting. 
There is some flexibility in choosing this function as long as y —■> 0 implies Ri —► 0. Some of 


the choices available may be understood using the kinematic variables described in Sec. 10 


If the secondary shower encoded in the Monte Carlo interface function to be described 


in Sec. |8.5| is based on the PYTHIA algorithm it might be sensible to set Ri to a scaled 
virtuality, 

Ri({p,f}m,y,z) = siy , (8.7) 

where si is the virtuality scale associated with parton l as discussed at Eq. (|6.4|). In the 
HERWIG case one might choose a scaled squared transverse momentum, 


Ri({p,f}m,y,z) = siyz(l-z) . 
We have indicated the scale of a s as a function 1st, defined by 

kr({p , f} m , l, y, z) = [si yz{ 1 - z)) l/2 . 


( 8 . 8 ) 


(8.9) 


Other choices for the scale in a s are possible. The matching to fixed order perturbation 
theory will work as long as we can write a s (kT({p , /} m , l, y, z)) = a s (pn)+0(a 2 ). Although 
our notation allows the possibility of a variety of choices for the functions Ri and fcr, for 
purposed of minimizing the dependence on the arbitrary parameter d\ n i our favored choices 
are given by Eqs. (|8.8| ) and (|8.9| ) with si a certain function of {p, f} m that is defined later 


in Eq. dl3T6|) . 

In Eq. ( |3.6D , we also make use of a function d defined in Eq. ( |6.4| ) that gives an 
approximate version of the resolution variable associated with the splitting generated with 
the splitting parameters yi , zi. With the use of this function, we limit the splitting in Ei 
to be unresolvable at a scale d that is approximately d m \ x 2 pi ■ Pkjsi (see Eq. Q6.5| )). 

In each Ei k operator, there is an operator on the parton color space, Ti ■ T k /[—T i J 2 ], 
that is familiar from the dipole splitting formulas. Next, there is an operator on the parton 
spin space, Si, with the interpretation that (a s /(27r)) Si is the probability for the parton 
to split at a given evolution parameter r if it has not split at a higher evolution parameter. 
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The splitting function Si depends on the splitting parameters Y] = {yi, Zi,(f>i, for 

parton l as well as on the momentum pi, which is needed to fully specify the meaning of 


4>i. We will specify this function later in Sec. 11, Eq. (11.4). 

The next factor, the Sudakov exponential, gives the probability that none of the par- 
tons has split at a higher evolution scale. Thus we work in a scheme similar to that of 
Sjostrand and Skands [14| and of Nason [15], picking out the hardest splitting. In the 
exponent there is a sum over partons l 1 and an integration over virtualities y' and the 
momentum fractions z' of the “virtual” splittings. The corresponding evolution parameter 
r' is required to be bigger than r. 

The remaining factor in the Sudakov exponent is the average over angle and flavors of 
S for parton 


{S{yv, = J 


dj>i 

2 7T 


x 

fl'l’fl 1 ,2 


j f f ~ 

fl'l+fl' 


Sv ( pv , yv, zi> , 4>i' ,fv, i , fv, 2 ) 


( 8 . 10 ) 


With our definitions, this is a numerical function times a unit operator on the partonic 
spin space. Explicit expressions are given later in Sec. |ll], Eq. (ED, (ED- 


8.5 Monte Carlo interface function 


The last factor in Eq. (8.1) 


is 




( 8 . 11 ) 


We imagine that after the parton splitting represented by the splitting function Ej which 
we may call the primary splitting, there is further parton showering, which we may call 
secondary showering. This showering is to be carried out by a shower Monte Carlo style 
computer program. The factor I represents the average value of the observable correspond¬ 
ing to the daughter hadrons after secondary showering when the shower starts with initial 
conditions specified by the variables ({p, /} m ; l, k, Yi). Here the initial conditions include 
not only the partonic state 


{p, f } m +1 = Ki <k ({p, f} mi Y{) 


( 8 . 12 ) 


generated from the partonic state {p, f} m by splitting parton l according to the splitting 
variables Y), but also the history of the first step of showering as specified by Y). 

We can think of I as being an integral, 


'=EpEn 

N {f}N i=1 


N 


dpi 5(p1 - m 2 (fi )) 
xP({P, f}N\{p, f}m ; l, k, Y t ) F N ({p} N 


(8.13) 


There is a sum over the number N of final state hadrons that are generated by the shower, a 
symmetry factor 1 /N\, a sum over the flavors {/}n of these hadrons and an integration over 
their momenta {p}n- In the next line there is a factor P that represents the probability 
density to produce the final state hadrons given the starting conditions represented by 
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({p, /} m ; l, k, Y[). The final factor in Eq. ( 8.13 ) is the measurement function evaluated 
with the produced final state hadrons. If one were to substitute F = 1, one would get the 
normalization condition for the conditional probability P , namely 1 = 1. Equivalently, if 
we substitute F = A, then 

| f=a = A . (8.14) 

Here A could depend on the initial conditions for the secondary shower as represented by 
the variables ({p, /} m ; l, k, Y t ). 

A simple model for / is obtained by omitting all secondary showering. Then I becomes 


I {0 X{pJ} m ;l,k,Y l ) = F m+1 (K l , k ({p,f} m ,Y l )) . 


(8.15) 


We will have more to say about the construction of the Monte Carlo showering program 


represented by function I in Sec. 15. One important feature is that the primary splitting, 
the splitting of parton l, should be the hardest splitting in the shower (according to the 
hardness measure R in Eq. ©)■ However, we leave the choice of this program largely 
open. Here we simply note two properties that / should have when the measurement F is 
an infrared safe N -jet measurement function. 

First, secondary showering should provide perturbative and power suppressed correc¬ 
tions to the simple function (assuming, always, an infrared safe observable): 


I{{p,f}m,l,k,Yi) = F m+1 (7l ltk ({p, f} m ,Y t )) 

x [1 + 0(a s ) + 0(1 GeV/Vs)] 


(8.16) 


Here the order cc s correction corresponds to splitting with a substantial virtuality, while 
the power suppressed correction corresponds to hadronization. 

Second, when yi vanishes, I should reduce to I ^ with only power corrections, 

I({p,f}m-J,k,0,z h (f>iJi,iJi,2) = F m ({p} m )x[l + O{lGeV/^)} • (8.17) 

The requirement here is that having the initial splitting virtuality yi equal to zero should 
set the maximum hardness for all of the secondary splittings to zero and thus turn the sec¬ 
ondary showering off except for hadronization. The hadronization model should then turn 
the partons into jets of a limited mass, leading to only power suppressed contributions to 
the measurement function. Note that F m+ i({p, /} m +i) reduces to F m ({p} m ) here because 
of the infrared safety property of the measurement function. 

In Secs. ^/7| and |3.8| , we will also need a function I{{p, f} n ) that represents the average 
value of the observable corresponding to the daughter hadrons after secondary showering 
when the shower starts with a partonic state {p, f} n with no other information (other 
than the d\ n [ cut) given as to previous shower history. The function I({p,f} n ) obeys the 
normalization condition 

I({pJ}n)\ F=x = \ , (8.18) 

where A could depend on the initial conditions for the secondary shower as represented by 
the variables ({p, f} n )- 
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As in the case of the version of / above, we assume that I is constructed so that the 
effects on the measurement of the showering are suppressed either by a power of a s or a 
power of 1 GeV/ y/s, 

/({p, f}n) = F n ({p } n ) x [1 + O(o s ) + 0(1 GeV/v'i)] . (8.19) 


8.6 Perturbative expansion 

We now seek the perturbative expansion of <r^ +s in Eq. (0- Insert into Eq. (|8.1| ) a factor 
1 = 7) + Ai, where 7] sets yi to zero in / and A i = 1 — 7). Then 

0 -m +S = ^ E / dr({p}m ) 0(d ini < d m {{p , f} m )) Wm({p, f}m ) 


{fh 


xEE (M({p,f} m )\ dY t E^ttpJU^lMdpJU)) 

1=1 k+l d 


( 8 . 20 ) 


x (Tt + A^IdpJU-J,^) . 


Our strategy will be to separate the 7} term from the A/ term. We call the the 7} term 


a 


and the A/ term Thus we write 


t B+S _ B+S 
T m ~ a m,T 


_B+S 

'm,A 


( 8 . 21 ) 


Then we will expand each term in powers of a s , up to next-to-leading order. We will 
find that, to this order, there are two terms in the expansion of 


a. 


B+S 

m,T 


B+S 


B+S 


Here crj^+ro} is proportional to a^ m and is, in fact, the Born contribution, , to a m . 
Then E+t+} consists of certain af 771+1 corrections. We will then find that the expansion 
of begins at order af 771+1 , so that we need only the first term in this expansion in 

order to evaluate <r5; +s to order a^ m+1 . 


a m,{T,0} + °m,{T,l} 


+ 0(af™ +2 ) 


( 8 . 22 ) 


8.6.1 The 7) contribution 

We begin by analyzing the term (J Y1+ ^. We use 

r|%/}m,i,fe,l|)) = %/}m;i,fc,0,«l,^ I /i ll ,/l,2) • (8.23) 

Using the property ( |S.17D of / we have 

H{pJ}m-,l,k,0,zi,<j ) i,fi, 1 Ji, 2 ) = F({p} m ) x [1 + 0(1 GeV/ y/s )] . (8.24) 

Then 

a m^T = —y E / dr ({p}m)0(d ini < dm({p, f}m )) W m {{p, f} m ) F m ({p } m ) 

m ' {fU J 

x EE (M({p,f}m) | fdYi EtMpJU^lMUpJU)) (8 - 25) 

Z=1 fc+Z J 

x [1 + 0(1 GeV/x/s)] • 
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Let us look at the sum of the E[ k operators integrated over the corresponding splitting 
variables. We have 



'' ^ rri rri 

\ ^ \ J-1 ■ J-k 
2^2^ _T' 2 In 

1=1 k^l l 


dr f — [ dziS(r - Ri({p,f} m ,yi,zi)) 
Jo m Jo 

a s (k T ({p,f}m,l,y,z )) 


X d(d({p, /} mi li Vh z l) < ^ini J 

x - V 8% - [ — Si{ Pl , fi,Y t ) 

2 fi,i+fi,2j 0 2ir iyyhJh u 


2t r 


fl,l,fl,2 


expf- [ dr' [ [ dz'^2 S(r'- Ri({p,f} m ,y',z')) 

\ Jr J 0 V J 0 // 

Oi s {k T {{p, f}m,V,yJ,z')) 


l> 

. ./ / \ 


X 9(d({p,f} m ,l,y ,z) <d h 
x {S(y',z',fi'))^ , 


2vr 


(8.26) 


where (S) was defined in Eq. ( |S.10 ). 
factors is 


According to Eq. (|4.6|), the sum over k of the color 


T l ■ T k 
Z. -j _T 2 

k ± l 


= 1 


(8.27) 


when operating on the state \M{{p,f} m ))- Having eliminated the color factor, we have 
an integral of a derivative, 



Rv({p,f}m,y',z')) 


X 0(d({p, f} m , l', y', z') < d\ n \) 
x (S(y',z',fi>))\ 


Us(kT({p,f}m,l',y\z')) 

2vr 


= 1 , 


(8.28) 


since the exponent is 0 at r = oo and — oo at r = 0 (since the integration region for r = 0 
includes the y —> 0 singularity). Thus (when operating on the state | M({p, /} m))) 



(8.29) 


We consider next the factor W m ({p, f} m ) in Eq. ( S.25| ) , which we write as 1 plus an 
order a s contribution cc s Wm {{p, f}m) according to Eq. (|jh3|), 


w m ({p , fU) = 1 + w£\{p, fU ) + ■ ■ ■ 


(8.30) 
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Taking the lowest order term in gives the lowest order term in the perturbative 

expansion of u^ +s . This term is 


a 


B + { T,°} = —j E / dr ({p}™) 9 { d ™ < d m({p,f}m )) 


X f}m)\M{{p, f}m)) Fm{{p}m) 

x [1 + 0(1 GeV/v/i)] ■ 

Comparing to Eq. (|6.2D , we see that 

<dT,o} = ^x [l + 0(lGeVA/i)] . 


(8.31) 


(8.32) 


The only order a* term in af n tj? comes from the expansion of the reweighting function 
W. This term is 


Sjr.i} = ~i E f dr ({p}m) d {dini < d m ({p, f} m )) F m ({p } m ) 

' ^^W^iipJj^iMiipJU^MiipJU)) 


{fir, 


x [1 + 0(1 GeV/Vi)] , 

where Wm is the first order contribution to W, Eq. (|8.3|). 


8.6.2 The A^ contribution 


Now we turn to the contribution in 
have 


Eq. ( |8.21 ). From the definition in Eq. ( 8.21 ) 


we 


E [ dr ({p}m) 8 {dim < d m ({p,f}m )) W m ({ P , f} m ) 

171 ' {f}J 

*EE (M({p,f} m ) | f dYi Ei tk ({p, f} m ,Yi) \M({p,f} m )) 

1=1 k^l d 

xa attpju-j,^) . 

Here 


AiI({p,f} m -J,k,Yi) = I({p,f} m -l,k,Yi ) - I({p, /} m ; l, k, 0, zi, <j>u fi,iji, 2 ) • (8.35) 


For / with yi = 0 we can use Eq. (8.24). For I with nonzero yi we use Eq. (8.16). This 
gives 


d{\Pi f}m'ih k,Yi) — |-F m +1 (F-lAiPJUrYt)) - F m {{p} m )} 

x [1 + 0(a s ) + O (l GeV/v^s)] • 


(8.36) 
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(8.37) 


We use this result in 


Eq. (|8.34ft , giving 


Ea = I dr ({p}m ) 0(d ini < dm({pj}m )) W m ({p, f} m ) 

m ' {fU J 

m 

(M({p, f} m ) | 

1=1 k^l 

X ^ Fm+ 1 

x [\ + 0{a s ) + 0(l GeV/y^)] . 


JdYt E ltk ({p, f} m ,Yi) \M({pJ} m )) 


As we will see, the leading contribution to crj^^ is of order af rn+l . Therefore, we 
can expand all of the factors and keep only the leading order terms. Using Eq. (|8.3| ) we 
can replace W by 1. We then use Eq. ( |8.6| ) for Eiy-, replacing the Sudakov exponential 
by 1, since we want only the first perturbative contribution (and since, because of the 
subtraction at yi = 0, the integrand is not divergent at y\ = 0). We can also replace the 
running coupling a s (kr ) in Ei^ by a s (pR) at leading perturbative order. Then 


= “j E / dr ({p}m ) d {dini < d m{{pJ}m )) 

' {f}m J 

m » 

AE dYi 1 7 Vh Zl) ^ ^ini) 

1=1 k^l J 

x <M({p,/} m )|^#Si(B,/i.r,)|x({p,/} ra )) 

X ^F m+1 {n l)k ({p, f} m ,Yi)) - F m ({p} m ) | 

x [l + 0{a s ) + 0(l GeV/ >/s)] . 


(8.38) 


In this equation there is an integration f dYi , which we can write out in full using 
Eq. ( 3.10 ). Then we can change integration variables to momenta {p } m +1 after the splitting 
described by Y). The jacobian is given in Sec. |I(J Eq. ( 10.17 ). After the change of variables, 
we have 


2- (m+ i)! E J dr ({p}m+l) 

V 7 {f}m +1 

xEE d ( d ({P’fYm k ’{hy’ Z }ij,k) < dim < dm({p, /}£’*)) 
k^i,j 

pairs 

x «s(m) 167T 2 

2tt 2pi ■ Pj (1 - y ijtk ) 

X (Map, /Jg*)> 

x{ F m+1 {{p}m + i)-F m {{p}^ k )) 

x [l + 0(a s ) + 0(l GeV/x/i)] • 


(8.39) 
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(8.40) 


Using Eq. 

T^ij,k ( {P ■ f }m+1) = 


— (M({p, m k ) \ y, z, <?%*, fi, fj)\M ({p, f }% k )) , 


2 Pi"Pi ' " ~T i} 

we see that if we identify 


fV 1 7-2 

Tr 1 -;- Si(p,f,Y) = Vi({p,y,z,4>},f 1 ,f 2 ) , 

2ir 1 — y 


(8.41) 


we will have 


B+S 


(m + 


l-j e f d nm m+1 ) E E 2\/,fc({P, /}m+l) 

{/}m+i hi Mhf 

pairs 


X 0(d{{p,fy^ k ,{l,y,z} ij>k ) < d ini < d m {{p,f }% k )) 
x {F m+1 (M m+1 )-F m ({p}^ fc )} 
x [l + 0(a B ) + 0(l GeV/v/5)] ■ 


(8.42) 


We will use in this form to combine with <7^ +s , which is discussed in the following 


section. 


8.7 NLO real emission corrections with shower 

We turn to the discussion of the NLO corrections. Let us start with the real contribution. 
Define 


= 

'- y 


1 


£ Jdrmm + l)l{{pJ}m + l)W^ +S ({p,f} m+1 ) 

{/}m+1 


X (M({pJ} m +i)\M({p, f} m +i)) 

X 6 ({j3, /}m+l) ^ ^ini ^ /}m+l)) 

-EE 'Dij.k ({p> f}m+ 1) 


(8.43) 


hi fc^hi 

pairs 


x 8(d({pJ}% ,{l,y,z}i j:k ) < d ini < d m ({p,f}% k )) 


This formula is illustrated in Fig. |j. The first term is the m + 1-parton matrix element 
squared and the second term is the sum of the dipole contributions to eliminate the infrared 
singularities. There is a Monte Carlo interface function / with the property Eq. ( 8.19| ). 
There is a reweighting factor W^; +S ({p, f} m +i ) that is similar to the reweighting factor W 
in cr® +S - We discuss this factor in Sec. [l3], Eq. ( |1 3.1 9| ). All that we need to know here is that 
W^ s ({p,f} m+ i) has a perturbative expansion that begins with 1: W^ +S ({p, f} m + 1 ) = 
1 + 0(a s ). There is a cut on d m and d m+ \. The m + 1-parton state should be not 
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resolvable at scale dim but, having put the two closest partons together, the resulting m- 
parton state should be resolvable. In the dipole terms we have cuts on d m and on d rn+ \ for 
the splitting in question, Eq. (|6.4|). The m-parton states defined by the dipole momentum 
set {p}m k should be resolvable at the scale d- m \. However, every dipole splitting gives only 
non-resolvable radiations according to the d measure. 


In Eq. (8.43) there is a potential infrared singularity when two of the partons, pi 
and pj, become collinear or one vanishes. This singularity is rendered harmless by the 
corresponding {ij, k} subtraction terms. Each of the {ij, k} subtraction terms can have 
other singularities when two of the momenta {p}m k are collinear or one is zero. However 


at such an extra singularity, d m ({p, f}]^ k ) vanishes, so that the cut dim < d m ({p, /}m A ) 
eliminates the singularity. 


1 ij,k\ 




Figure 2: Illustration of c^ +s , Eq. ( |3.43 ). As in Fig. |l], the three partons that cross the dash lines 
are resolvable at scale dj n i. In the first graph, there is an additional perturbative splitting that is not 
resolvable at scale di n i. There is a subtraction graph in which the perturbative splitting is replaced 
by a splitting according to the dipole splitting formula, which is the order a s contribution to the 
splitting represented by the box in Fig. [l]. The subtraction removes the soft/collinear divergence 
from the splitting. Finally, there is secondary showering represented in the figure by diamonds 
and in Eq. ( 8.43| ) by the function /. There is a reweighting factor W that is not represented by 
a graphical symbol. This illustration is an attempt to depict the main idea behind Eq. ( |S.43| ). 
Actually, though, the full matrix element \M) to produce four partons appears in the first diagram 
and the full matrix element \M.) to produce three partons appears in the second diagram along 
with a sum over the dipole subtraction terms. 


The perturbative expansion of this is simple. Using Eq. ( 8.19|) , we can replace / by 
the jet observable function F and, using Wj^ +S ({p, /} m +i) = 1 + (D(a s ), we can replace W 
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by 1. This gives 


(T R+S = 

w m 


1 


(m + 1 )! 


dr ({p}m+ 1) F m +1 {{p}m+ 1) 


{/}m+l 

X f} m+ i)\M{{p, f}m+l)) 

X $ (^m+1 ({]3, /}m+l) ^ ^ini ^ /}m+l)) 

-EE Aj,fc({P,/}m+l) 

ij fc+J 
pairs 

x 6(d({p,f}% k ,{l,y,z}i jjk ) < d ini < d m ({p, f}% k )) 
x [l + 0(a B ) + 0(l GeV/v/i)] . 

Now if we add ct r+S to from the previous section, Eq. ( 8.42| ), we get 


R+S, B+S 
a m + (T mA 


(dr+ E /^({#) m+ 0 

{/}m+l 

(•Ad({p> /} m+ i)|A4({p, /} m+ i)> ^m+l({p}m+l) 

X 9(d m +i ({p, /}m+l) r?i n ; < dm({P) /}m-t-l)) 

-EE A j , t ({p,/} m+ i)F m ({p}«’ t ) 

i,j k^i,j 
pairs 

x 0(d({p,/}^{Z,z/,4v,k) < d ini < d m ({p,/}^’ fe )) 
x[l + 0{a s ) + 0(1 GeV / i/s) 

Comparing with Eq. ( |iO|) , we see that 

1 + 0(a s ) + 0(1 GeV / \/s) 

8.8 NLO virtual corrections with shower 

In this subsection, we turn to the virtual loop corrections. We define 5 


_R+S I _B+S _ R-A v 
V ®m, A X 


u v+s = 
w m 


J2 [ dr ({PU) V{ d ini < d m ({p, f } m )) I({ P , f} m ) W^ +S ({ P , f} m ) 
m ' {fU J 




l k^l 


^^W^H{pJ} m )(M({pJ} m )\M({p,f} m ))\ . 


5 The expression in Eq. (8.47[) can be simplified using Eq. (|l2.8|). 


(8.44) 


(8.45) 


(8.46) 


(8.47) 
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This formula is illustrated in Fig. ||. We integrate over the phase space for m partons with 
a cut on d m ({p. /} m ). There is a Monte Carlo interface function I with the property ( 8.19| ). 
There is a reweighting factor W^ +S ({p, /} m ), which we discuss in Sec. [L3], Eq. ( 13.19[) . This 
function has a perturbative expansion that begins with 1: W^ +S ({p, f} m ) = 1 + 0(a s ). 
Then there is a factor with three terms. The first two contain the functions V and C) 
from the perturbative virtual loop contribution, Eq. (^T^). The third contains the first order 
contribution Wm^ to the Sudakov factor and is similar in structure to Eq. ( 8.35 ). 



Figure 3: Illustration of cr^ +s , Eq. (8.47). As in Fig. |l], the three partons that cross the dash lines 
are resolvable at scale dmi- To the left of the first dashed line, we illustrate one contribution to the 
one loop amplitude | M^). The infrared poles of | M^) are subtracted according to the dipole 
subtraction scheme. To the right, we illustrate one contribution to (A4|. All of this is included in 
Eq. ( |g.47 ) in the function V, which is described in Sec. 12. Finally, there is secondary showering 
represented in the figure by diamonds and in Eq. (8.47) by the function I. There is a reweighting 
factor W that is not represented by a graphical symbol and there are two additional terms in 
Eq. (8.47|) that are not illustrated in the figure. 


To find the perturbative expansion of this, we replace I by the jet observable function 
F, using the property (|S.19|) of I and, using W^ +S ({p, f} m ) = 1 + 0(a s ), we replace W 
by 1. This gives 


= — 

^ m i 

ml 


y. / dr({p} m ) 9(d mi < d m ({p,f } m )) F m ({p} r 


{/T 


V({p,f}m) ~ EE Cl, k ({ P }m,d ini ) 


l k^l 


(8.48) 


- f}m) (M({p , f}m)\M({p, f}m )> 

1 + 0(a s ) + O (l GeV/\/s) 
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We now add cr® 4 ^ 1 | from Eq. (|3.33|) 


V+S I q-B+S _ 
a m w CT mi { T ,l} — 


^2 [dr({p}m ) 6»(d ini < d m ({p,f}m )) F m ({p } m ) 

EE Cl,k({p}m,d in i) 

l k^l 

1 + O(os) + 0(1 GeV/ y/s) 


(8.49) 


Comparing with Eq. ( |Q|) , we see that 

= °l +A [l + 0(«s) + 0(1 GeV / y/s) ] . (8.50) 

8.9 Net result 

Combining the contributions cr® 4 ^ Q | from Eq. (|8.32| ) , cr^ +s + from Eq. ( |8.46| ) , and 
v+s + 212\T 1} f rom Eq. ( 8.50| ) we see that the calculation of a m with showers, using both 
the Born term and the correction terms, reproduces a m to next-to-leading order accuracy, 

(8.51) 


o, 


B+S , R+S , V+S B , R—A , V+A 
i '-'m ' 


+ 0(af m+i ) + 0(a* m x 1 GeV /y/s ) 


This is the result claimed in Eq. (7.2). 

We also note that the calculation of a m with showers using only the Born term repro¬ 
duces a m to leading order accuracy, 


<7m +S = 2a + 0(af m+1 ) + e>(af m x 1 GeV/Vs) . 


(8.52) 


9. Alternative LO partial cross sections with showers 


In Eq. B we imagined that the perturbative partial cross sections a m are known at 
order af m+1 but that beyond a certain value ttinlo of rri only order af m results are known. 
Furthermore, we suppose that exact order a,f m results are known only up to certain value 
TO mav of m. We needed a construction of the corresponding contributions cr^ 4S including 
showers. Our choice was to simply define crj)) +s f° r m > m NLO according to Eq. ©■ This 
reproduces cr m to leading order. 

In this section, we wish to point out an alternative. For m > tunlo one could replace 
o"m +S by ct^ +s defined by 

d^ +& = | dr({p} m ) 6(d ini < d m ({p, f } m )) W m {{p, f}m ) 

(9-1) 

X f}m)\M({p, f}m)) i({p,f}m ) ■ 


Here we use the m-parton Born cross section modified by the reweighting factor W with a 
cut di n i < d m ({p , f} m - We omit the dipole style splitting that we used in Eq. (8.1). Instead, 
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we pass the partonic state {p, f} m to the shower Monte Carlo program, represented here 
by the Monte Carlo interface function I with the property Eq. (|8.1S), as used already in 


Eqs. (8.43) and (8.47). The shower Monte Carlo program generates parton showers from 
the scale di n i downwards. 

This is to say that for m > tunlo we simply use the scheme of Ref. [11]. This is 


perhaps a little simpler than the use of a splitting based on the dipole subtraction scheme 
for the first splitting. 

There is a largest value, m max , of m such that we have available exact leading order 
matrix elements. Thus the sum in Eq. ( 0 ) stops with cr^ +s or <r® +s for m = m max . This 
leaves us with a zero cross section for producing more than m max jets that are resolvable at 
scale G?i n i. Even though we lack the matrix elements for calculating exactly at leading order 
the cross section for producing more than m max jets, we can still have an approximate cross 
section for producing many resolvable jets. Assuming that we are using cr® + , all that we 
need to do is make a simple modification in the formula fl9.1|) for <t^ +s for m = m max : 
in W m ({p, f} m ) and in I({p, f} m ) we should replace d ini by d m ({p, f} m ). That is, there 
are splittings that produce m max partons according to the exact matrix element, with the 
corresponding Sudakov factors and a s factors, represented in the modified W m ({p, f} m ). 
We remove from W m ({p, f } m ) the factor representing the probability that there was no 
splitting between d m ({p, f} m ) and d; n i because we now allow such splittings as part of the 
modified I({p, /} m ), which describes an indefinite number of splittings at scales below the 
scale d m ({p, f} m ) of this last splitting according to the exact matrix element. 


10. The kinematics of parton splitting 

In this section, we review the dipole splitting construction of Catani and Seymour 
using the notation adopted for this paper. This covers one parton splitting into two partons 
(with the participation of a spectator parton). With a trivial extension that we have already 
described in Sec. |j, this construction also covers one of m partons splitting to create an 
m + 1-parton state. 

The construction and deconstruction of dipole splitting outlined in Sec. § is based on 
a transformation 

{pi,Pj,Pk} <->■ {Pij,Pk,y,z,4>} • (io.i) 

The transformation from left to right combines partons i and j. With the inclusion of the 
rather trivial transformation to combine the flavors, this transformation was called Qij^k 


in Eq. (T8). The transformation from right to left splits parton ij and, with the inclusion 
of the flavor splitting, was called IZij (or IZi) in Eq. (|3.7| ). 

Combining parton momenta works as follows. Starting with three massless momenta 
{PiiPjiPk}, Catani and Seymour define pij and pk by 


Pij =Pi+ Pj - 
Pk = Pk + 


Pi-Pj 


(.Pi +Pj)'Pk 
Pi-Pj 


Pk , 


( 10 . 2 ) 


(Pi+Pj)-Pk 


Pk 
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Then it is evident that 

and that 


Pij +Pk=Pi+ Pj + Pk 


Pij =Pk = 0 • 


The splitting parameter y is defined by 


V = 


Pi-Pj 


Pi-Pj + (Pi +Pj)-Pk 


(10.3) 

(10.4) 

(10.5) 


so that 


Pij = Pi + Pj ~ 


i - y 


Pk 


l 


( 10 . 6 ) 


Pk = 


i - y 

The splitting parameter z is defined by 


Pk 


z = 


Pi-Pk 


(Pi+Pj)-Pk 


(10.7) 


Then (1 — z) is given by the same expression as for z but with i j. Alternative formulas 
for y and z are 

Pi-Pj Pi-Pk nne x 

y = - —v- , ^ = -—— • (io.8) 

Pij'Pk Pij'Pk 

To define the azimuthal angle (j) we first define k± to be the part of (fp — pj )/2 that is 
orthogonal to both p v} and Pk, 


k± = \ ( Pi - pj) - (z - \)pij + {z-\) ypk 


(10.9) 


One can show that k± obeys k± -p^ = 0 and k±-pk = 0 using the easily proved relations 
(pi — Pj)-pk = (2 z — l)pij-Pk and (f) t — Pj)-Pij = (2 z — 1) ypij-pk■ A convenient alternative 
formula for k± is 

fei = zpi - (1 - z)pj - (2z — l)pij . (10.10) 

The squared length of fcj_ is 


k]_ = —2y z(l - z) pij ‘Pk 


( 10 . 11 ) 


The unit vector 



( 10 . 12 ) 


defines the plane of the splitting in a reference frame in which the vector parts of pij and 
pk lie along the positive and negative z-axis, respectively. The azimuthal angle of k± with 
respect to some convenient reference direction is (j) ■ Stated perhaps more precisely, we 
say that one specifies 0 as a shorthand for saying that one specifies the unit vector k±. 
Integrating over 0 means integrating over n± subject to the conditions that it is a unit 
vector orthogonal to p VJ and pk- 


- 34 - 



The inverse transformation, giving {pi,Pj,p k } in terms of {pij,Pk,U, z,<j>}, is easily 
obtained by combining Eq. ( 10.9 ) and Eq. ( 10. 6|) , 


Pi = zpij + y(l - z)p k + [2yz{l - z)pij-p k \ 1/2 
Pj = { 1 - z)Pij + yzpk -[2yz{l- z)p ir p k \ vl2 K ± 
Pk = { i - y)pk ■ 


(10.13) 

It is reasonably straightforward to work out the jacobian for this transformation. Defin¬ 


ing 


dr (n \{p} n , Q) = f[ ((27 r)- 3 d 4 Pl 8+(pt)j (27r) 4 5 4 - Q^j , (10.14) 


one finds 112 ] 


dr^(pi,pj,p k -,Q) = dr^\piyp k ;Q) dydz ^ (1 - y) 


2tt 16tt 2 
x 0(z(l - z) > 0) 0(y(l - y) > 0) , 


(10.15) 


or, equivalently, 


167T 2 


d,r (i> (pi, p r /p k : Q) 1U/ ‘, = dr^ipijiPk-Q) dydz ^ -—- 
2p*-Pj 2 tt y 

x 0(^(1 - z) > 0) 0(y(l - y) > 0) 


(10.16) 


This transformation between two parton momenta plus splitting variables and three 
parton momenta is trivially extended to a map between m parton momenta plus splitting 
variables and m + 1 parton momenta, as described in Sec. ||. If the parton that splits is 
labelled l, the jacobian is 

dr({p} m ) dyi -——dzi = dr({p} m+1 ) , . (10.17) 

yi 2vr 2pi t i ■ p i>2 

In our applications, this result will appear with symmetry factors and sums over the 
indices of splitting partons. Then it looks like 


m 


^ Y J dr ({p}m ) EE 

‘ {/} m J i=i kjtl 

rl d yi f 1 ^ 


xi — dzi 
I o Vi Jo Jo 


$yE tf u+fi2 (i-v,)H{pj} m+1 ) (1018) 




1 


(m + 1) 


; Y / dr ({p}m+l) Y Y 


{f}m +1 


hi k+ij 

pairs 


167T 2 
2 Pi ■ Pj 


H{P,f}m+ 1) , 


where h({p, f} m +i) represents any well behaved function. To prove this, note that the 
left hand side has symmetry factors 1/m! and 1/2 and a sum over m values of the index 
l and (m — 1) values of the index k. On the right hand side we have a symmetry factor 
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l/(m + 1)! and a sum over (m + l)m/2 values of the index pair i,j and (m — 1) values of 
the index k. If we put the sums over parton indices outside the integrals, each term is the 
same, so we can just take one term and multiply by the number of terms. We note that 
the net counting factor m(m — l)/[2m!] on the left hand side equals the net counting factor 
(m + 1 )m(m — l)/[2(m + 1)!] on the right hand side. Thus Eq. ( |10.17| ) implies Eq. (10.15). 

We now examine the relation between parton splitting and the jet resolution functions 


dij defined in Eq. 


and d defined in Eq. 


Taking massless parton momenta in 


Eq. (541) and using Eq. ( 10.5 ), we have 


, „ x 2 pi- pj . 

d^ = d(pi,pj) = -— mm 


pi-n pj-n 
pj-n pi-n 


2pij-p k . I Pi-n pj-n 

y mm< —-, —— 

Pj-n pi-n 


(10.19) 


where n is a unit vector that defines the time axis in the e + e _ c.m. frame. For pi-n and 
Pj-n, we can use Eq. ( |10.13|) . Of particular interest is the case 0 < y <C z < 1, in which pi 
is either collinear to pij or both soft and collinear. Under this condition we have 


r 2 pij -pk 

dij ~ —-- y mm 


1 - z’ 



2 Pij-Pk j 

s ij 


( 10 . 20 ) 


Here d is the function defined Eq. (|6.4D (with yi —> y, zi —> z and with the scale si for 
parton l renamed to s^). We used this function to limit parton splitting in Eq. (|6.3p and 
Eq. (|8.6|): d < d\ n \. We see that, under these conditions, limiting d limits dij , although 
there is a factor 2pij-pk/sij that relates the size of the two resolution measures. 

Also if interest is the case 0 < z < 1, 0 < y < 1, i//z ~ 1, in which pt is soft but not 
collinear to pij or pk- Under this condition, we have 


, 2 pij-p k 

dij ~ —-- y z 


l + V Pk-n 
z pij-n 



( 10 . 21 ) 


Thus in the soft limit, dij has the same scaling behavior as d, i.e. a factor yz, but the ratio 
of these functions is not generally 1 and depends on the ratio of y to z. 

We should also be concerned about the region 0 < z <C y < 1. This region is not en¬ 
hanced in the integration over y and z because the relevant splitting functions in Eqs. ( |6.3| ) 

0 at fixed y. Nevertheless, since d -C dij in this region, 
means. Consider, then, the resolution 


and (S.6) are not singular for z 
we should understand what the restriction d <C d m \ 
parameter for combining parton i with the spectator parton k. Using the definition (| 
together with Eq. ( |10.8|) , we have 


d ik = d{pi,p k ) = 


Oij-Pk x . f Pi-n Pk-n 

^ -v I I n i \ m i vj^ / 


z (1 — y) min 


Under the condition 0 < z <S y < 1, this is 

, 2 Pij-Pk . L (i -y ) 2 

d ik - yzmmi 1, -^— 

s { y z 


{Pk-n Pi-n 
2pij 'Pk 


< 


d . 


( 10 . 22 ) 


(10.23) 




Thus when 2 <C y the condition d < di n i restricts parton i to be close to the spectator 
parton k. 
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11. Splitting functions 


In this section we review the functions used in the dipole subtraction algorithm |]T2| . 
We then specify the functions Si(p,f,Y ') and ( S(y,z,f )) that we use to describe parton 
splitting and that are derived from the T>ij^ functions. 

We begin with the functions 'D%j.k that form the basis for the dipole subtraction algo¬ 
rithm. The basic idea is that the squared matrix element for 

m + 1 partons has a potential singularity that could lead to a divergent integral when the 
dot product of any pair of the momenta, say pt ■pj , goes to zero. The matrix element takes 
a rather simple form in this limit. In fact, the form would be extremely simple were it not 
for the fact that one has pi ■ pj —> 0 not only when pi becomes collinear with pj but also 
when pi—> 0 with pj fixed and pi not necessarily collinear with pj (or when pj becomes soft 
with fixed pi). This configuration can lead to a divergent integral if the soft parton i (or j ) 
is a gluon. When the soft gluon i couples to parton j and another parton k, the structure 
of the limiting function depends on what k is. For this reason, Catani and Seymour write 
the matrix element in the limit pi ■ pj —> 0 as a sum of terms labelled by the index k of a 
“spectator” parton. That is 


(M({pJ}m+i)\M({pJ} m+1 )) = ^ V ijik ({p, f} m + i)H-, (11-1) 

k^i,j 

where the dots stand for terms that are nonsingular in the limit pi-pj —> 0. The dipole 
functions 'Dij k have a simple structure of the form that we reviewed in Eq. ([f.5|), 


f}m+ 1) — 

1 


2 Pi ■ pj 


<vi({p,/}S‘)| ry T ‘ 


_ rp2 


ViA\pu-n- fin*, n, fj) | M({p, m k )) 


( 11 . 2 ) 


First, there is a singular factor 1 /(2 pi ■pj ). Then there is the Born amplitude | M({p, f}m k )) 
and its complex conjugate for rn parton momenta and flavors formed by combining partons 
i and j with the help of the spectator k according to the formulas of the previous section. 
These amplitudes are vectors in the color and spin space of the partons. More precisely, 
for each of the m partons there is a spin space spanned by two basis vectors | s) and a color 
space spanned by three basis vectors |c) in the case of a quark or antiquark or eight basis 
vectors |c) in the case of a gluon. The amplitude \M.({p, f} m )) lies in the direct product 
of the m spin spaces and m color spaces. In Eq. ( 11. 2| ) there are color and spin operators 
that act on given single parton factors in the direct product space, with the parton factor 
affected labelled by a subscript ij (for the parton obtained by combining partons i and j ) 
or k (for the spectator parton). 

For the color, there is an operator, the SU(3) generator T^“, that acts on the color 
space for parton ij and there is another SU(3) generator Tjj that acts on the color space 
for the spectator parton k. The dot product indicates a sum over a from 1 to 8. In the 
denominator there is a factor Tjj ■ T^, which is simply a number equal to Ca if parton ij 
is a gluon and Cj? if it is a quark. 
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The remaining factor, Vj, is a function of variables defined by considering that one 
combines partons i and j with the help of the spectator k, as in Eqs. 


|) and (fPp. 

Specifically, Vij is an operator on the spin space of the mother parton before the splitting 
and depends on the momentum pij of this parton. It is also a function of the splitting 
variables {y, z,cf>} and the daughter flavors 

The shower algorithm of this paper makes use of splitting functions Si(p, Y ). We found 
it useful to take these functions to be proportional to the functions Vij , using Eq. (|S.41), 


a s 167T 2 
2tt 1 — y 


Si(p,f,Y) = V l ({p,y,z,(j)}J 1 j 2 ) , 


(11.3) 


where Y = {y, z, 4>, /i, ^ 2 }- (Note that the subscript l does not enter into the functional 
dependence of these functions but merely tells on what parton’s spin space the operator 
acts.) Since we use mostly the functions S[(p, /, Y), we present here the standard definition 
of V.j simply translated into the new notation. 

The definition begins by separating the possibilities for flavors, 


Sl{p,f,Y) = S gfl Sgf 2 Sgg(j),Z,y,K ±) + ^ 

r=u,u,d,d ,... L 


5 r fi $rf2 Sqq {p, Z,y, Kj_ ) 


+ drfa 6gf 2 Sqg ( Z,y ) + dgfr 5 r f 2 Sqg (1 z, y) 


(11.4) 


Then for the splitting of a quark or antiquark into the same flavor quark or antiquark plus 
a gluon the splitting function is 


{s I Sgg(z,y) | s') = Cf (1 - y) 


_1 - z( 1 - y) 


-(1 + *) 


(11.5) 


where the s and s' are the spin indices of the emitter quark or antiquark. For the splitting 
of a gluon into a quark and an antiquark, we again denote the spin indices of the emitter 
gluon by s and s' and define 


(s\S q q(p,z,y,K±)\s') = Tr (1 — y) e^p^) [-g^ - Az(l - z) k^_k u ± \ e u (p,s') 
Finally for the splitting of a gluon into two gluons the splitting function is 
(s\Sgg(p,z,y,K±)\s') = 


( 11 . 6 ) 


2C a (1 — y) el(p,s) - 


+ 


l-z(l-y) l-(l-z)(l-y) 
+ 2z(l - z) k^/c^ e v {p,s') . 


- 2 


(11.7) 


In the gluon splitting functions, e(p, s ) is the polarization vector for the emitter gluon. With 
a change of the gauge used in defining e(p, s ), one has e(p, s ) —► e(p, s) + Xp. However, the 
matrix elements are unchanged because e(p, s) ■ p = 0 and k± ■ p = 0. 

It is also of interest to know whether these matrix elements depend on the spectator 


momentum, p No spectator momentum appears explicitly in Eqs. (11.6) and (11.7), but 
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recall that k± is a unit vector orthogonal to p and to pf~. If we change pi- to p' k , we can 
define a new vector k [|_ to specify the azimuthal angle of the splitting by 


/ Pk ' K -L 

K ± = -^- p 

Pk ' P 


( 11 . 8 ) 


This new vector is still a unit vector, still orthogonal to p, but now is orthogonal to p' k 
instead of pk■ The change does not affect the matrix element because e(p, s) ■ p = 0. This 
justifies the notation that Vij depends on p and on the splitting variables (including </> or, 
equivalently, k±) but not on the spectator momentum pk- 

The Sudakov exponent in the showering formula contains the average over angle and 
flavors of Si, Eq. Q8.10 ), 


(S(y, z, /)> =f^ ■ ( 1L9 > 

flj2 


Here the angular average is in 2 transverse dimensions. There are extra terms in ( S ) if one 
works in 2 — 2e dimensions, but we do not need these terms in this paper. For the angular 
average, we can use 


/ 


dcj) 


^(p,s) K ± K ± e v(P ’ s 0 



( 11 . 10 ) 


Then we see that {Si(y, z, /)) is a number times the unit operator on the partonic spin 
space. A simple calculation gives for a quark or antiquark emitter 


( S(y,z,f)) = C F (1 - y) 


.1 - z(i - y ) 


-(i + *) 


/e{u,u,d,...}. (li.il) 


For a gluon emitter, including splittings into both gluon and quark-antiquark pairs, one 
finds 


< S(y,z,g )) = C A (1 - y) 


+ 


l-z(l-y) 1 — (1 — z){l — y) 
+ lVfT R (1 - y) [l-2.z(l-;z)] . 


-2 + 2(1 - 2 ) 


( 11 . 12 ) 


The y —■> 0 limits of these functions match the Altarelli-Parisi splitting functions, 

( S(y,Z,f )) -»• C F = P q/q ( z ) : / G {u, U, d, . . . } , 

(S(y, z, g)) C A [1 ~ < Z 1 (1 _~, Z)]2 + N { T R [1 - 2z(l - 2 )] (11-13) 

2(+ Z) 

= \ P g/g( z ) + N t P a/g( z ) ■ 

Note that some of the Altarelli-Parisi functions are singular as 2 -> 0 or 2 -> 1. However, 
for y > 0 the splitting functions S are not singular as 2 —> 0,1. 
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12. Virtual contributions 


In Eq. (4.7), we specified the virtual loop contribution in the dipole subtraction scheme in 
terms of a function V({p, /} m ), 


a 


V+A _ 


[dcr Y + [ da A } = -1 £ [ dr ({p}m) V({pJ}m ) F m ({p} m ) . ( 12 . 1 ) 

Jl m ' {f}m J 


The definition [12] 


is 


V({pJ} m ) 


2 Re ( M (M m )| M {1) (Mmi e)) 

+ Re (M({p} m )\l({p, f}m'i e)|M(Mm)) 


( 12 . 2 ) 


Here ({p} m ;e) ) denotes the virtual one-loop correction to the m-parton matrix el¬ 

ement in 4 — 2e dimensions. This matrix element is singular as e —► 0. The singularity 
is cancelled by the integral of da A , performed in 4 — 2e dimensions. The result of this 
integration is the Born amplitude | -Ad (MM) times an operator 1(e) in the color space of 
the final-state partons, 




Cfco 


27t T(1 — e 


-i m m m m / a 9 

1 ^ ■ T.j ( 47 Tfi 


2=1 


(12.3) 


j¥* 


The function V(f, e ) has a simple expansion about e = 0, 

V(/,<0 = c, (4 - y) + 7/ (1 + l) + Kf + 0(c) . (12.4) 

where Cf, 7 /, and Kf are given in Eq. s 

When we adapt the dipole subtraction scheme to the calculation of partial cross sec¬ 
tions based on a cut c?j n j on the “distance” between partons, we need a correction to the 
subtraction term for virtual graphs. The correction involves a function C) ^ defined in 
Eq. Q, 


Ci,k({p,f}m,dmi) = j ^ J dz J ^^Y^S f ^ + ^d(d- m i<d({p,f} m ,l,y,z)) 


/1J2 

w y(i - v) 2 Pr Pk ^ , r ^ ^ 

x - ^2 - P l,fc({pj}m;^) • 


(12.5) 


Here the distance function d({p, z) is defined in Eq. (|6.4| ). To find the explicit 

form of Ci t k({p, f} m , dmi)j we use Eqs. ( 11.2 ), ( 11.3 ) and ( |11.9| ) to relate the sum over 
daughter flavors and average over azimuthal angle 4> of T>ij. to (S(y, z, /)). This gives 


Ci, k ({p,f}m,d ini ) = gc( (M({p,f} m )\^^\M({p,f} m )) , (12.6) 
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where 


«s 

27T 


r 1 dv r 1/2 

C{x,f) = — dz0(z/(l-z) >x/y)[(S(y,z,f)) + (S(y,l-z,f))\ . (12.7) 

Jo y Jo 


When we sum over the spectator label k, Eq. ( 12.6| ) simplifies to 

s d[ n [ 


Y C l,k{{P , f}m, d-mi) = ^ C 

k^l 


Si 


,fl (M({p, f}m)\M({p, f} m )) 


( 12 . 8 ) 


Using the explicit form of (S(y, z, /)) given in Eqs. ( |11.11| ) and ( |11.12| ), we can perform 
the integrations over y and z to obtain 

1 t r 2 

C(x, f) = -C f log 2 (a;) + 7 / log(4x) - 2 C f log(x) log(l + x) - — C f 

+ 7 / 2 xlog( 2 .x) - 27^(1 +x)log(l + x) + 7/(1 - x) - 2Cf\A 2 {-x) (12.9) 

1 + 2x — 3x 2 ' 


/1 2 x 

+ S g f(C A - 2 T r N f ) ^-xlog —— + 


12(1 + x) 


13. Sudakov factors 


The functions u® +s , and it^ +s contain reweighting functions W. These functions 

contain rather trivial ratios of a s at different scales and not so trivial Sudakov exponentials. 
In addition, the splitting operators E/.k in <r® +s contain Sudakov exponentials. In this 
section, we discuss these factors and, in particular, give definitions for the functions W 
that we use. The presentation also touches on the Sudakov exponentials in the the Monte 
Carlo interface functions / and I in <r^ +s , <r^ +s , and <r^ +s . 


13.1 The standard NLL Sudakov factor and its interpretation 

The formula (8.1) for <r^ +s contains a reweighting factor W m ({p, f} m ) consisting of a 
product of ratios of o; s at different scales, Eq. (8.2), and a Sudakov exponential. Following 
Ref. |llj, we base this factor on a construction of a splitting history corresponding to 
the final state { p , f} m using the “hr” jet algorithm (slightly extended to include flavor 
information) as described in Sec. j5| The algorithm combines the final state partons so 
as to produce a QCD tree graph for e + e~ —► {p, f} m . Let the lines in this tree graph 
be labeled by an index L. We can label the vertices according to the order in which the 
partons were combined. We start with m partons and combine two of them at a vertex 
with label V = m to produce an m — 1 parton state. Then we combine two of these m — 1 
partons at a vertex with label V = m — 1 to produce anm-2 parton state. Finally, the 
last 2 partons are combined at vertex 2 to form a 7 or Z. At each stage between vertices 
V and V + 1 there are n = V partons. The line labels L for these partons form a set that 
we can call I(V) (for the intermediate state after vertex V). 

Each vertex V has an associated distance measure d(V), with d(V) > d(V + 1). The 
values of the d(V) will appear in W m ({p, /} m ). Typically, d(V) is the hr measure (|5.3|) for 
the splitting at vertex V. In some cases, we have redefined d(V) to equal d(V + 1). 


- 41 - 














With these definitions, for each line L in the splitting graph there is initial resolution 
parameter d(Vj(L)) equal to d(V) for the vertex Vj (L) at which the line originates and a final 
resolution parameter d(V{(L)) equal to d(V) for the vertex d(Vf(L)) at which the line splits. 
For a line that enters the final state, we take d(V{(L)) = di n i. Thus 1 > dfL) > df (L) > d- in i. 
We denote the flavor of line L by f(L). 

With this notation, we can state the definition of the reweighting factor W as given in 
pd| 1. Briefly, the idea is the following. Consider the Born cross section er® with observation 
function F m ({p} m ) = 1. This is the Born level cross section for making m jets, using 
resolution parameter d; n i■ By inserting a factor W m ({p, f} m ) into the integral for <x®, we 
obtain the approximation to the ?n-jet cross section that sums logarithms of di n ; at the 
leading log and next-to-leading log level. The required reweighting factor is a product of 
a s factors, one for each vertex, and a product of Sudakov factors, one for each line, 


< LL (o>,/}„) = n a ‘ { VWTs) a m m{L)),<k, 

- LJ - ™ (Mr) 


V >2 


a. 


A f(L) (d(Vf(L)),d ini ) 


(13.1) 


where 


Af (d m ax> d m i n ) — exp 


f dmax dX a s (V\s ) f / dmax 

L. T^^w log (— 


7/ 


(13.2) 


The coefficients Cf and 7 / are given in Eq. ©• Thus one multiplies the perturbative 
matrix element by Sudakov factors that would have been associated with the lines if the 
final state had been generated by a suitable parton shower algorithm according to the 
specified QCD tree graph with the restriction that no splittings were to be generated 
that were unresolvable at resolution parameter d[ n [. For each QCD vertex we also use the 
reweighting factor to replace o; s (/xr) used in the perturbative matrix element by the running 
coupling used in a standard parton shower. Our default reweighting factor W m ({p, f } m ) 
will be a slightly modified version of W^ LL ({p, f} m ) that is equivalent if we neglect terms 
suppressed by a color factor 1 /N%. We define the modified version below after examining 
the structure of VF,^ LL ({p, f} m )- 

The function W^ LL ({p, f} m ) as given by Eq. (|13.1 ) is a product of factors, one for 
each propagator in the jet-structure graph. There is another way to write it that focuses 
on the evolution from one value of d(V) to the next (See Fig. [I]). 

To reorganize B / )^ LL ({p, f} m ), we define a special vertex V S (L) corresponding to each 
line L in the graph. We can think of V S (L) as the vertex where the history of line L 
started (cf. Sec. 3.2 of |Tl|). The identification of the V S (L) is recursive. Consider the 
starting vertex 2 : 7 jZ —> q(Li)q(L 2 ), where our notation indicates that 2 is the label 
of the vertex, L\ is the label of the quark line and L 2 is the label of the antiquark line. 
We define V S (L\) = V S (L 2 ) = 2. Then for any vertex V) : q(Lo) —► q(L\)g(L 2 ), we define 
V S (L) for the two daughter lines by V S (L\) = V(Lq) while V S (L 2 ) = Vi. Similarly for any 
vertex V) : q(Lo ) —> q(Li)g(L 2 ), we define V S (L\) = V s (Lq) while V S (L 2 ) = V). For a vertex 
Vi : g(Lo) —> g(Li)g(L 2 ), we first provide a definite labelling for the daughters by, say, 
letting L\ be the label of the daughter gluon with splitting fraction z\ > 1/2. Then we 


- 42 - 












3 



Figure 4: Illustration of the tree graph for a shower constructed from a five parton final state. 
The lines in the graph are numbered L = 1,... 8. The vertices are numbered V = 2,..., 5. The 
resolution parameters d(V) for the vertices decrease: d( 2) > d( 3) > d( 4) > d( 5) > dmi- There is a 
starting vertex V S (L) for each parton line L. Here 14(1) = 14(2) = 14(3) = 2, 14(4) = 14(5) = 3, 
14(6) = 14(7) = 4, and 14(8) = 5. We will write the Sudakov factor as a product of factors that 
represent the probabilities that the partons in each of the indicated intermediate state evolve from 
scale d(V) to scale d(V + 1) without splitting. For the last final state, the Sudakov factor is the 
probability for the five partons to evolve from scale d(5) to d; n ; without splitting. 


define V S {L\) = V{Lq) while 14 (^ 2 ) = 14- Finally, for any vertex V) : g(Lo) —> q(Li)q(L, 2 ), 
we define V S (L\) = V s (L 2 ) = V S (L 0 ). 

In the case that L labels a quark or antiquark line, we define another special vertex 
V 0 (L) that is the vertex where the quark line originated in the jet-structure graph: either 
the 7 /Z —* qq vertex V = 2 or else a later g —> qq vertex V. (The subscript is “o” for 
“originated.”) 

The reader can check that with this definition 


W* LL ({pJ} m ) = n Yl exp{-S nll (V,L)} 

. (mr) 


V>2 


(13.3) 


V=2 LeI(V) 


with 

^NLL 


f 

(V,L)= / 
J d. 


d(y) dA a s (\ // As) / „ . fd(V s (L))\ 


d(v+1) , 


A 27t 1 

+ 6»(/(L) G {«,IZ, d,... }) 


— -C F 
2 


logfffiMV 

S \d(Vo(L))J 


(13.4) 


Here we interpret the lower limit d(V + 1) to be di n ; in the case V = m. This represents 
a simple algebraic reshuffling of the original expression that appears in Eq. ( |1 3.1 ). What 
is more significant than the algebra is the interpretation. The total Sudakov factor is a 
product of factors, one for each splitting vertex V. Each of these factors represents the 
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propagation of the system between scales d(V) and d(V + 1). The factor for propagation 
from d(V) and d(V + 1) is a product of factors, one for each parton line L that occurs 
in the intermediate state after splitting vertex V. Each of these factors, exp{— S{V, L)}, 
represents the probability that parton L in the intermediate state following vertex V did 
not split at a resolution parameter between d(V) and d(V + 1). In the case of the last 
vertex, V = m, exp{—<S(V, L)} represents the probability that parton L in the intermediate 
state following vertex m did not split at a resolution parameter between d(m) and d m \. 

In the Sudakov factor, we integrate over the resolution parameter A = k^/s of a virtual 
splitting. As we will see, we can understand the logarithms in the integrand as the result 
of integrating over the angle of a virtual splitting: 


. fd(V s (L))\ [d(V B {L))s/Ei dQ 2 

) = L/e; 

/ d(V s (L)) A _ rd(V4 L ))s/E* d0 2 
l0g \d(V 0 (L))J ~ 

Jd(V 0 (L))s/E2 02 

, f d(V s (L))\ fd(V s (L))\_ fd{V 0 (L))s/ES dd 2 

g ( A ) ,0g {d(V 0 (L))J J„ /Ei <P ' 


(13.5) 


Here E s is the energy of the least energetic of the two partons that were produced at vertex 

V a (L). 

The angle limits that appear in these formulas can be understood on the basis of 
“angular ordering” |L6[]. There are three types of splitting vertices to consider, g —a gg , 
q a qg , and g —a qq. Then the result in Eq. (|1 3.4 ) can be understood as the result of 
stringing together the three cases in all possible combinations. 

The case of a g —> gg vertex, V, is easiest. Typically, one of the daughter gluons is 
“hard” and moves in almost the direction of the mother gluon. The other gluon is softer 
and moves at an angle By = kxy/Ey with respect to the mother direction, where Ey is 
the energy of the softer daughter gluon. That is, 


0 2 y = d(V)s/Ey. 


(13.6) 


Suppose that we now emit a soft gluon with energy E from the two daughter partons. The 
soft gluon has transverse momentum kx with respect to the daughter parton from vertex 
V. The resolution parameter for the soft gluon parameter is A = kj,/s and we integrate 
over A in Eq. (13.4). The emission angle 6 of this soft gluon is given by 6 2 = A s/E 2 . The 
smallest that 6 2 could be is (9j( lin = A s/E 2 mx . where we can take E max to be the energy 
of whichever parton is the source of the emission. Now suppose that the emission angle is 
less than the angle By of the g —> gg splitting. Then we get incoherent emission, with color 
factor Ca, from each of the daughters. If we emit a soft gluon at a larger angle from the 
two daughters, the gluon does not resolve the daughters and instead sees the color of the 
mother. Thus large angle emission comes with a color factor Ca- Thus we get soft gluon 
emission with a color factor Ca for either small or large angles with respect to the hard 
daughter direction, and we get additional radiation with a color factor Ca for small angles 
with respect to the direction of the soft daughter gluon. 
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For emission from the soft daughter, with a factor Ca, we have an angular range 
d(V)s/Ey > 0 2 > A s/Ey, where we have set E max = Ey. In this case, vertex V is the 
starting vertex for the soft daughter, so that the angular range is d(s)s/E 2 > 9 2 > A s/E 2 , 


as in the first line of Eq. (13.5). 

For emission from the hard daughter plus coherent emission from both daughters, we 
have, to start with, a maximum angle equal to the angle of the vertex at which the mother 
gluon was a daughter gluon in a splitting. However, if the mother was the hard daughter 
parton in a g —► gg splitting, we can also have coherent emission involving the mother’s 
sister. Continuing in this way, we get emission with color factor Ca up to a maximum 
squared angle 9 2 where vertex s is the vertex at which the mother of the hard daughter 
started as a (relatively) soft daughter of energy E s (either in a g —> gg splitting or in a 
q —> qg splitting). That is 9 2 = d(s)s/E 2 . The minimum angle is given by 9 2 rim = A s / E 2 xax 
where we can approximate the hard daughter energy E max by E s . Thus the angular range 
for emissions from the hard daughter is d(s)s/E 2 > 9 2 > A s/E 2 , as in the first line of 
Eq. fll3lip . 

What about a q —> qg vertex? Typically, the daughter quark is “hard” and moves in 
almost the direction of the mother quark. The daughter gluon is softer and moves at an 
angle ky/E with respect to the mother direction. Suppose that we emit a soft gluon from 
the daughter partons with an emission angle less than the angle of the q —> qg splitting. 
Then we get incoherent emission, with color factors Cy or Ca respectively, from the quark 
and gluon daughters. If we emit a soft gluon at a larger angle from the two daughters, 
the gluon does not resolve the daughters and instead sees the color of the mother. Thus 
large angle emission comes with a color factor Cy. That is, we get soft gluon emission with 
a color factor Cy for either small or large angles, and we get additional radiation with a 
color factor Ca for small angles. 

These two cases account for the factor log(d(F s (L))/A) that appears in the first term in 
Eq. (13.4). To understand the remaining term in Eq. (|13.4| ), we need to consider a. g qq 
vertex. 

In the case of a g —> qq vertex, both of the daughter partons are typically quite hard. 
We have incoherent soft gluon emission with color Cy from the two daughter partons as 
long as the emission angle is less than the splitting angle at the g —> qq vertex. For larger 
angles, we have soft gluon emission with color factor Ca because the soft gluon sees the net 
color of the qq pair, which is the color of the mother gluon. It is convenient to ascribe half 
of the large angle emission probability to each of the quark and the antiquark daughters. 
That gives an emission probability Ca/2 for each. We can understand the Sudakov factor 
in Eq. ( 13.4 ) in the case that f(L) £ { u , u,d,...} as describing an emission probability Cy 
for emissions from the quark line with angles smaller than the angle of the qq vertex at 
which the quark line was created and with a factor Ca for larger angles, up to the angle 
of the vertex V S (L) at which the gluon that split to make the qq pair had its start. 

13.2 Modified W for the Born cross section 


Now that we understand the physical picture behind Eqs. (13.3) and (13.4), we can propose 


a small modification. As was pointed out in [11] using rather different language, the 
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term in Eq. ( |13.4 ) proportional to [Ca/2 — Cf] may seem unnecessarily complicated and 
could perhaps be eliminated. Note that this contribution arises only when there is a 
g —*• qq vertex, whereas g qq splitting is rather unlikely, particularly because the g —> gg 
comes with two logarithms while the g qq splitting comes with only one logarithm. 
Furthermore, the color factor is small compared to Cf: [Ca/2—Cf]/Cf = 1/(N 2 — 1) = 1/8. 
Shower Monte Carlo programs typically drop contributions that are of relative order 1 / N 2 , 
so it seems appropriate to do so here. For these reasons, our default choice is to drop this 
term and instead define W m ({p, f} m ) using 


W m ({pJ} m )= [I Qs (v^>) exp {-S(V,L)}, (13.7) 

n (Mr) 


V>2 


V=2 Lei(V) 


with (using d(m + 1) = d;, 




Since we include the factor W in the perturbative matrix element in <r^ +s and we wish 
to maintain the NLO accuracy of the calculation, we need to keep track of the first term 
in the perturbative expansion of W, Eq. (|8.3|), 


w m ({ P , f } m ) = i + w£\{ P , f} m ) + 


Using Eqs. ( 13.7 ) and ( 13.8 ) and using Eq. 
we have 


(13.9) 

to express a s (y/dvs ) in terms of a s (/di), 


w£\{p,f} m )=-J2 E 

V =2 Lel(V)' 

m 

- E log 

y=3 


(W dX 

d{V+ 1) ^ 

d(V)s' 


C f(L) log 


d(Vs(L)) \ 


V(L) 


Mr 


(13.10) 


where, again, d(rn + 1) = dj n ;. 

13.3 Sudakov factors in dipole splitting 

Let us now look at the Sudakov exponential that occurs in the function Eq. (|3.G|). 
We set the evolution variable to be k^/s, as in Eq. (|8.8| ). Here = siyz( 1 — z) and kx 
is also the argument of the running coupling according to Eq. (^1]). Then the Sudakov 
exponential in Eq. (^) is independent of l and k and is 

A(r) = JJexp^— J dr 1 J J dz 1 5{r' — si'y'z!(1 — z')) 


x 6 


r 


max-jV 2 , (1 — z') 2 } 

a s {Vr i ) 


< dini'S 


(13.11) 


2vr 
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Here the product runs over all final state partons l'. After using Eqs. (|11.11|) and ( |11.12j ) 
for the splitting functions, we can perform the z'and y' integrals in the exponent to obtain 


A(r) 


n ex p 

V 



dX a s (\/Xs) 

A 2ir 


+ o 




+ 0 


C f(i') lo § 


si'/s 

A 





(13.12) 


Here sy/s lies between d m i and 1. This form for the exponent applies when A <C s^/s, with 
a correction that vanishes when A/(s^/s) —> 0. There is a second correction term that is 
bounded by a constant times 


h(d ini / A) = 6 (A > d ini /4) log(l/(2[l - v / V4i])) • (13.13) 


This correction comes from the difference between having a simple cut r'/s < dim and 
having our more complicated cut r'/s < dim x max-jV 2 , (1 — z') 2 }. The two correction 
terms make finite contributions to the Sudakov exponent when r/s -C sy/s with dim either 
small compared to sy/s or of order sy/s (but in any case bigger than r/s). 

We see that, to the accuracy of dropping non-logarithmic terms in the exponent, we 
have 


A(r) ~ exp 


f dini d\ a a (VXs) j ( si/s\ 

l T 2ir ] C,m l 0 g V ~) - 


(13.14) 


Compare this to the factor in Eq. ( 13.7 ), corresponding to the probability that the final 
state particles not split between the scale d(m) and the limiting scale dim. According to 
Eq. (13.8), this is 


n 


0 —S(m,l) _ 


n«p 


i i 

Let us choose 


^ m UX a s (VX^) 


X 2vr 


C f(D lQ g 




W) 


Si = S dy(£) 


Then 


A( r ) n e 


—S(m,l) _ 


n exp 


^dX a s (V As) \ 
X 2ir Y 


Cm . 0 Y-m) 


X J 


T7(0 


(13.15) 


(13.16) 


(13.17) 


corresponding to the final state particles not splitting between the scale d(m) and the scale 
r, at which one of them does split. There is no dependence on dj n j in this product. 

We thus see the purpose of including the factor W rn ({p. f} m ) in the formula for cr® +s : 
without it, there would be a dependence on d; n i from the Sudakov factors for showering 
below the resolution scale d; n i. More precisely, this is the purpose of including the factor in 
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W m ({p, f}m ) relating to propagation from scale d(m ) to dim. The other Sudakov factors in 
W m ({p, f}m ) are equally important. The complete expression for cr® +s has Sudakov factors 
from three sources. Evolution from c ?2 to d m is included in W m ({p, /} m ). Evolution from 
d m to r = d m+ 1 is included partly in W m ({p, f} m ) and partly in the Sudakov factor A(r) 
that is part of Eik, as we have just seen. Then there is a succession of shower splitting 
scales that we have represented as being part of I({p, /} m ; l, k, Yj). This showering will 
involve Sudakov factors for evolution from m + 1 to m + 2, from m + 2 to m + 3, and so 
forth. If we hold all of the parton momenta fixed and change dj n ; from just less than d m to 
just greater than d m . then the same event history counts not as part of <7® +s but rather as 
part of i • The integrand representing the probability for this event history in cr^ +s will 
be somewhat different from the integrand for the same event history in <7^^, but the two 
integrands will be approximately the same. The main points are that the Sudakov factors 
from all three sources are approximately the same and that the m + 1-parton squared 
matrix element is approximately the m-parton squared matrix element times the dipole 
splitting functions. 


13.4 Factors W for real and virtual corrections 


We have defined the factor W m ({p, f} m ) in <7® +s and examined its structure and how it 
meshes with the Sudakov factors in the splitting function Ei^. We now turn to the function 
^ / m +S ({P'/}m) ^at appears in cr)( t +s , Eq. ( |3.47 ), and the function W^f({p, f} m+ i) that 
appears in cr^ +s , Eq. ( |S.43|) . 

For IFj)( +s ({p, f }m) j we have an m-parton final state that is constructed to be resolv¬ 
able at scale dim. It thus seems sensible to use the same reweighting function that we used 
for the Born contribution, 


WZ +S ({pJ}m)=W m ({p,f} r 


(13.18) 


It is worth noting that this choice is not as obvious as it was in the for Recall that 

cr^ +s is based on graphs with a virtual loop (together with a subtraction that keeps the 
virtual partons from being very soft or collinear with the external particles). If the virtual 
partons always had momenta of order s/s then the choice ( 13.18[) would be physically well 
motivated. However, the momentum scales in the virtual loop can be intermediate between 
\fs and the scales of the parton splittings in the synthetic shower history constructed 
from the final state {p, f} m . For this reason, it could well be that the “best” choice for 


W^ +S ({p, f}m) is something more subtle than that given in Eq. ( 13.18 ). 

For Wm +S ({p, f}m+i), we have an m + 1-parton final state, so that we have resolution 
scales d( 1) > ■ ■ ■ > d(m ) > d(m + 1) in the synthetic shower graph generated from the 
final state {p, f} m + 1 - In the main term of er,^ + , we have d(m) > dim > d(m + 1). There 
are then subtraction terms that remove the leading singularity when d(m) 3> d(m + 1). 
For this reason, d(m + 1) is typically not much smaller than d m i, nor is d(m) much larger. 
We therefore define W^ +S ({p, f} m + 1 ) to include a reweighting factor for a s at each strong 
interaction vertex and a Sudakov factor giving the probability that the partons did not 
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radiate between each pair of vertices: 


m+1 


w^ +s ({ Pl f} m+ 1)=n 

V=3 


a s (\/d(V)s) 

■S ( mr ) 


n n ex p{- 5 ( F > L )} 

V=2 Lel(V) 


(13.19) 


with S(V,L ) as defined in Eq. (13.8). The last Sudakov factor takes us from vertex m, 
after which there are m partons, to vertex m + 1 at which one of these partons splits with 
a scale d(m + 1) to make an rri + 1-parton state. Subsequent evolution, as given in the 
Monte Carlo interface function I K+S ({p, f} m+ i) in Eq. (|S.43D , then starts at splitting scale 
d(m + 1). Again, it is worth noting that this choice is not as obvious as it was in the for 
<x R+s because of the subtraction term in <r R+s . It could well be that the “best” choice for 
Wm +S ({p, f}m) is something more subtle than that given in Eq. ( 13.19 ). 

An alternative is to take 


W 

’ ► 7 


■V+S/ 


({p,f}m) = W* +S ({p,f} m+ 1 ) = 1 


(13.20) 


This choice removes some sensible physics built into W^ +S and W R+S . However, it does 
have a technically useful feature. Suppose that we take a measurement function F = 1, so 
that cr R+s + <r^ +s + cr R+s is the contribution to the total cross section from final states 
with precisely m jets. The Born contribution, <r R+s , contains terms of all orders in ct s 
starting at order af 771 . It does not have the complete contribution at order a^ m+1 and 
beyond, but it does have the correct leading and next-to-leading logarithms of di n ; [|ll|] (at 
leading order in 1 /N%, since we have modified W m a little). The addition of < 7 ,)() +s and 
(j R+s fills in the missing pieces at order a f 771+1 . It also can add contributions at higher 
orders. However, if we choose W^ +S ({p, f}m) = P^m +S ({lh /}m+i) = 1> then cr^ +s and 
cr R+s are exactly proportional to af 771+1 , with no contributions at higher orders. They 
provide just what is needed at order a^ m+1 and no more. To see this, note that in the 
expression (|3.43 ) for ct r+s , the factor I is 1 when F = 1 because of the property ( 8.1 8] ) 
of I. When we also set VE R+S ({p,/} m ) to 1, we get an expression proportional to a Rm+1 
with no further dependence on a s . Similarly, the choice H / J ^ +S ({p, f} m ) = 1 makes 


Eq. ( 8.47 ), exactly proportional to a^ m+1 . 


14. What size is the resolution parameter? 

The algorithm described in this paper depends on a jet resolution parameter d in i, which 
plays a role similar to that of the factorization scale in calculations of cross sections for hard 
processes with one or two hadrons in the initial state. We anticipate that d in i would be a 
parameter in a computer code that implements this algorithm, so that a user could choose 
its value. In this section, we describe some considerations that would go into making the 
choice. 

Suppose that we calculate u[F] for an infrared safe IV-jet observable that measures 
large momentum scale features of events and hence does not involve large logarithms. An 
example for N = 3 is F = T 4 , the fourth moment of the thrust distribution, Eq. O). 
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More generally, we can consider the nth moment of the trust distribution for n > 2, 

c[T n } = . ( 14 . 1 ) 

Given that we want to calculate cr[F] for an observable in this class, we seek to choose the 
resolution scale dj n ; so the dominant contribution to cr[F] comes from a m [F], with m = N. 
Let us examine why one would like ctn[F] to be dominant and how dominant it should be. 

We first consider the comparison of ctn[F} to <7jv+i [F], supposing that a m [F) for m < 
N — 1 is negligible. We note that the cross section cr[F] = Ylm^rnlF] has a perturbative 
expansion that starts at order atf N and has higher order contributions. For a next-to- 
leading order calculation, we need the next contribution, of order af N+1 . The partial 
cross section ctn[F] also has a perturbative expansion that starts at order af N . As we 
have seen, the algorithm presented in this paper reproduces the order af N and a ^ N + 1 
terms of this expansion. The perturbative expansion for cqv+il-^] begins at order a^ N+1 . 
Thus part of the order af N+1 correction to cr[F\ is in <t^ +s [F] and <t^ +s [F] and part is in 

| [F]. (Part is in a m [F} for m < N — 1, but we suppose for the moment that this part 
is negligible.) 

Consistently with this counting of powers of a s , we expect <T 7 v + i[i ? ] ~ a s <7jv[-F]. How¬ 
ever, suppose that we were to choose such a small value for di n i that a s log 2 (l/di n i) is of 
order 1. Then the cross section <jn[F] that we hoped was dominant would be suppressed 
and the cross section would be shifted to cr m [F] with m = N +1, N + 2 ,.... This would not 
much matter at leading order since the for-jet matching scheme arranges that the sum of the 
<7 ,^ +s [F] is very accurately independent of d in i. However, it is possible that the perturbative 
correction terms of order a^ N+1 can be lost among terms like af N+1 x ot s log 2 (l/di n i). Such 
a term is beyond the next-to-leading perturbative order for <tjv[-F] and beyond the next- 
to-leading logarithm approximation also. Terms like this are always present in [ [F]. 
(j^r +s [F] and but they are numerically important only when d - m i is very small. To 

ensure that these contributions are not numerically important, we can ask that the relation 
between ctn[F] and cttv-i-i[T 71 ] follow the perturbative expectation, 

cbv+i[-F] < a s £Tat[F] . (14.2) 

With a computer implementation of the algorithm described here at hand, one can also 
check that cr^ +s [F] and er^ +s [F] are of the right perturbative size (that is, a s aiy[F]) and 
that the complete &[F] is accurately reproducing the purely perturbative result for the same 
quantity, which can be obtained by simply turning off the shower part of the program. 

Consider next the partial cross section cy/v-i [T], with N — 1 resolvable jets. The 
perturbative expansion for <7Ar_i[-F] has the form 0 

<tn-i[F\ = C N . lfi [F]af N - 1 (Q) + C N . 1A [F]a^(Q) + --- . ( 14 . 3 ) 

This compares to the expansion for aj\r[F], 

a N [F] = C n , 0 [F] a^ N (Q) + C N}1 [F } a* N+1 (Q) + • • • . ( 14 . 4 ) 
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The first perturbative coefficient for ow-if-E] vanishes since, by assumption, F(f) = 0 for 
a state with N — 1 partons. Fortunately, we have the next term , 6 proportional to af N . 
We do not have the complete coefficient proportional to af N+1 . Thus it seems that we do 
not have the af N+1 accuracy that we want. However, F(f) = 0 for a state with N — 1 
very narrow jets. Furthermore (Tn~i [-F] is the contribution from final states with N — 1 
jets that cannot be further resolved at scale d- m \. If we take d; n i to be small, the N — 1 jets 
are guaranteed to be narrow. We conclude that all of the coefficients Cn~ij[F] will be 
small if we choose to be small. Since we are missing one factor of a s in the expansion 
of ai\r-i[F], we can simply demand that d; n i be chosen so that the coefficients Cn~ij[F] 
are smaller than the corresponding coefficients Cn j- \ [F] in cjn[F] by a numerical factor 
that is approximately equal to a s or else smaller. That is, we ask that 


<jn-i[F ] < a s a N [F] 


(14.5) 


Can one really satisfy Eqs. ( |14.2| ) and ( |14.5|) ? Consider the three-jet observable T 4 , 


Eq. ( 14.1 ). Using Pythia with y/s = Mz to estimate cr m \F^\ we find that <73 [T 4 ] is largest 
near d; n 1 = 0.04. At this value of dj n ;, we have 04 [Z4] = 0.08 <73 [T 4 ] and 04 [T 4 ] = 0.12 03 [T 4 ], 
thus satisfying Eqs. ( 14.2 ) and ( 14.5 ). 

The case of the second moment of the thrust distribution, F = T 2 , is instructive. 
Again taking yfs = Mz , we find that 03 [T 2 ] is largest near d m i = 0 . 02 . At this value 
of dim, we have <74 [T 2 ] = 0.16 <73 [T 2 ] and 04 [T 2 ] = 0.35 03 [T 2 ]. In this case the two jet 
contribution is not as suppressed as one might like. We do not interpret this as indicating 
a difficulty with the h ?-jet matching scheme. Rather we take it as an indication that, after 
accounting for showering and hadronization, cfXjj] gets a 30% contribution from the two 
jet region (defined with jet resolution 0 . 02 ) even though T 2 is nominally an infrared-safe 
three-jet quantity. For this reason, a pure NLO fixed order calculation is not as reliable 
a calculation as one might like and one should not be concerned if the algorithm in this 
paper based on kx -jet matching does not exactly match the pure NLO calculation. To the 
extent that they turn out to differ, we expect that the kx ~jet matching calculation with 
some perturbative accuracy sacrificed in favor of showers and hadronization is the better 
calculation. 

These examples suggest that a choice of d- m i in the range < di n i < a s might be 
appropriate in practical cases. We can also ask whether a choice of d m \ in the range 
ctg < di n i < a s would work in the asymptotic limit a s —> 0 . 


Consider first the condition in Eq. ( |14.5| ) . It is pretty straightforward to see that for an 
infrared safe IV-jet observable F, < 7 at_i[.F] will vanish like a power of di n ; (times logarithms) 
for di n i —► 0. For example, for the nth moment of the thrust distribution, Eq. & , one 
finds a limiting form 02 [T n ] oc d”{ 2 . Thus if we take di n ; ~ a s or smaller, Eq. ( 14.5 ) will 
be satisfied asymptotically as long as n > 2. Following the philosophy stated above for 
02 [T 2 ] when a s = a s (Mz), we might take the observables T n for n < 2 to be sufficiently 
contaminated by two jet physics that we would not need to satisfy Eq. for those 

observables. 


6 Since the required NLO calculations are available for N jets, we can assume that they are available for 
N — 1 jets. 
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Consider next the condition in Eq. ( 14.2| ). Since the perturbative expansion of (Tat-|_i[-F] 
begins with one more power of a s than the perturbative expansion of cr^vf-E], Eq. ( 14.2| ) 
appears to follow. There is a potential problem from the Sudakov exponentials, which 
could suppress <7/v[E] relative to er/v_|_i[T]. The exponents contain factors a s log 2 (l/di n i). 

a s or d- m i ~ a 2 , the Sudakov exponent becomes const. a s log 2 (l/a s ). In the 


0, these factors get smaller. Thus the Sudakov exponentials do not change the 


With dini 
limit a s — 

perturbative power counting and Eq. (14.2) is satisfied. 

Thus both numerical examples based on Pythia at a s = a s (Mz) and an analysis of 
the limit a s —> 0 suggest that a choice of d; n i in the range a 2 < d; n ; < a s might be useful. 
When one has a working program, one will want to revisit these questions with a numerical 
analysis using the program. 


15. The Monte Carlo interface function 


In Sec. 3J3, we introduced a function I({p, /} m ; l, k, Y{) that represents the expectation 
value of the observable in hadronic states that are generated by a shower Monte Carlo 
event generator beginning from initial conditions represented by the partonic state {p, f} m 
and the information about the splitting of parton l contained in the variables l, k, Y\. This 
function is used for <r® +s . We also introduced a similar function I({p, f} n ) that is used for 
ct^ +s and cq^ +s (and may optionally be used for <r® +s for m > tunlo)- For the matching 
to a NLO calculation, we needed only certain basic properties of these functions, given in 
Eqs. ( |S.16| ) and Q8.17D , and for /, in Eq. ( |S.19| ). 

In this section, we delve a little further into the requirements for the shower Monte 
Carlo event generator. There are quite a number of successful shower algorithms available 
and it is beyond our scope to specify any algorithm in detail. However it may be useful to 
say a little about the initial conditions for the shower Monte Carlo program. In particular, 
we note that we are handing the Monte Carlo a partially developed shower and that the 
Monte Carlo needs to begin at the point where the previous shower simulation left off. 

Consider, then, the function /({p, /} m ; l, k, Y\) used for er^ +s . There is an m-parton 
state that can be thought of as having been created by a showering process with a hardness 
cutoff d ini < d m ({p, f} m . After that, one of the partons (labelled l) splits. Suppose for the 
sake of concreteness that this dipole splitting is based on the evolution variable of Eq. (EH1), 


r = si yz( 1 - z) 


(15.1) 


Let us define a transverse momentum k± for this splitting according to Eq. (10.9). With 
this choice, we have 


\k±\ = 2pi ■Pk yz{ 1 -z) = 


2pi • Pk 
si 


(15.2) 


Note that for a given value of the evolution variable r, \k'^_\ for the splitting of parton l 
with the participation of spectator parton k depends on l and k. 

Now, all of the partons are allowed to split, and the one that does is parton l. The 
others did not split at an evolution variable above the value r. That is, parton V, with the 
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aid of spectator k', did not split with r' > r. Further evolution of these partons should 
then be restricted to the range r' < r. That is 


\k'l\ < 


2pv ' Pk' s l 


si> 2pi ■ p k 


\kl\ 


(15.3) 


A restriction like this can be imposed in the chosen shower Monte Carlo program by using 


a veto algorithm, as described for instance in Ref. [111. The only problem may be that the 
Monte Carlo program may know the index l 1 of the parton that it is proposing to split, but 
may not have a choice for a corresponding spectator parton k!. A sensible choice would be 
to let k' be one of the final state partons to which parton l' is color connected (at leading 
order in 1 /N c ) according to a rule based on the synthetic splitting diagram obtained by 
applying the kx- jet algorithm to the m-parton state and adding the one splitting of parton 
l. (For a gluon V there are two such color connected partons and one would choose either 
of them with probability 1/2.) An alternative that avoids selecting a spectator parton is 
to impose 


I A; | I < max 
k' 


2 pv ■ Pk' si 


\kl\ 


si> 2pi ■ p k 

For the splitting of one of the daughters of parton l, one may simply impose 


(15.4) 


\k'l\ < \ki\ 


(15.5) 


Of course, if the splitting angle is not the evolution variable used in the shower Monte 
Carlo program used for secondary showers, then the program should separately arrange for 
angular ordering. 

Consider next the function I({p, f} m ) as used for o// +s . The final state, which is the 
initial state for the Monte Carlo program, consists of m partons that can be thought of 
as having been created by a showering process, with a hardness cutoff d\ n [ < d m ({p. f} m . 
There is no further splitting before this state is passed to the Monte Carlo program. There¬ 
fore further splittings l —► i + j should be generated with a cut d(pi,pj ) < d\ n [. The same 
reasoning applies in the case that one elects to use / without a dipole splitting for cr^ +s i n 
the case m > uinlOj as described in Sec. |9| 

The case of I({p, f} m +i) as used for <r^ +s is similar. Here there are m + 1 partons 
with d m+ i({p, f} m + 1 ) < dini■ Thus it is d m+ i({p, f} m + 1 ) rather than d m \ that should serve 
as the upper cutoff for splittings included in I({p, f} m + 1 )- 


16. Conclusion 

We have proposed an algorithm for adding showers to next-to-leading order calculations 
for e + + e _ —» N jets. This algorithm incorporates two kinds of matching, which we might 
call PS/NLO matching and m-jet/(m + l)-jet matching. 

For the m-jet/(m + l)-jet matching, we make use of the kx -jet matching scheme of 
Catani, Krauss, Kuhn, and Webber |Ll|. We divide the cross section <y[F] corresponding 
to an observable F into contributions a m [F ] corresponding to final states with m hard jets. 
Here “hard” is defined using a jet resolution parameter dj n In the case that F is an IV-jet 
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observable that is infrared-safe in the technical sense and is in fact insensitive to physics 
at momentum scales much below y/s, we want the dominant contribution to a[F] to come 
from a m [F] with m = N. 

For PS/NLO matching, we are concerned with the m partons that emerge from a 
Born graph that contributes to a m [F]. One of these partons can split as the first step of a 
parton shower based on this graph. On the other hand, there are perturbative corrections 
to the same graph that also involve real or virtual splittings of the same m partons. The 
MC/NLO matching is needed in order to incorporate the parton shower and the NLO 
corrections without double counting. With the parton shower in place, we not only get 
<j m [F\ correct to NLO, but we also get sensible answers to more detailed questions about 
the low momentum scale structure of the m jets that are used to calculate a m [F]. A pure 
NLO calculation cannot (sensibly) answer such questions. 

We thus have a good description of hard physics and soft physics and even both at 
once. What about physics at an intermediate momentum scale? Suppose that F is such 
that a[F] is sensitive to physics at a scale K 2 that is similar to d; n iS- An example would 
be the thrust distribution da/dt with K 2 = (1 — t)s similar to dmiS- In this case neither 
( 72 [F] nor <73 [F] dominates the calculation of cr[F]. If we hold K 2 fixed and increase 
d; n i a little, <t 2 [F] decreases and <73 [F] decreases. If we use the Born level contributions 
<7® +s [F], the approximations used in i7^ +S [F] and cr^ +s [F] are so similar that the sum 


<7. 


B+Sr 


[F] + <73 [F] is very insensitive to d; n j. If we now include cr^ +s [F] + < 7 ,^ +s [F], we 

improve the approximations by adding NLO correction terms. However, the correction 
terms [T 71 ] + u^ +s [F] and ojf +s [I ? ] + cr^ +S [F] are not matched in the same sense 
that the leading terms < 7 ^ +S [F] and < 7 f +S [F] are matched. To do better at matching the 
correction terms for different numbers of jets would require next-to-leading order splitting 
kernels for the showers, which we do not have. For this reason, one may expect that adding 
the correction terms cr^ +s [F] + a^ +s [F] introduces some d m \ dependence into the result. 
Fortunately, the NLO correction terms are small, so the added di n i dependence must also 
be small. 

We have been quite specific about the subtraction scheme to use in the base NLO 
calculation. We have been deliberately less specific about the scheme for showering. The 
hardest splitting of a parton from a Born graph is to be done using the dipole splitting 
functions. Softer splittings can be handled by the reader’s favorite shower Monte Carlo 
method as long as the initial conditions for this secondary showering are compatible with 
the organization of the part of the showering already completed. 

It may be useful to compare the algorithm discussed in this paper with those in other 
papers. We can compare treatments of the Sudakov factor Ei^ in Eq. (8.1), the reweighting 
factor W, and the factor I that represents further showering. 

Consider first the exponent in the Sudakov factor F/^ for the splitting that is to be 
matched to NLO perturbation theory, Eq. ( | 8 . 6 | ). In the work of Frixione, Nason, and 
Webber js], f|, this exponent is the splitting function of HERWIG. Thus one must subtract 
something from the NLO calculation that is derived from the order a s expansion of a 
HERWIG shower. As discussed by Nason (15], this means that the NLO calculation is to 
be tuned to the shower Monte Carlo. Additionally, since the HERWIG shower does not 
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precisely match the perturbative probability for soft wide-angle gluon emission, one must 
also provide a special treatment for the soft region. The algorithm of Ref |8], ||] is based 
on the simple idea that, in a physical gauge, the cut two point function represents parton 
splitting, while the virtual two point function represents the probability not to split. For 
this reason, the Sudakov exponent is chosen as the integrand of the virtual two point 
function in the Coulomb gauge. There is a separate Sudakov exponential for emission of 
soft wide-angle gluons. This method works, but is not easily adapted to the dipole style 
of NLO calculations that most people use. In the present work, as in the leading order 
program of Sjostrand and Skands 0^ the hardest splitting in a shower gets a special 
treatment. Nason [jljj has proposed a scheme of this style that is adapted to an angle 
ordered shower, as in HERWIG. In this scheme, the Sudakov exponent for the hardest 
splitting is the ratio of the exact matrix element (A4({p, f} m +i)\AA({p, f} m + 1 )) for m + 1 
partons to the matrix element {A4({p. f} m )\AA({p, f} m )) for m partons, before one of the 
partons is imagined to have split. The m + 1-particle matrix element is to be divided into 
a sum of terms associated with individual parton splittings using an algorithm that is not 
completely specified in (l^] . There is a certain conceptual simplicity to this choice, but it 
puts a very complicated function into the Sudakov exponent. In this paper, the Sudakov 
exponent is the angular average of the splitting function in the dipole subtraction scheme 
that is widely used for NLO calculations. With this choice, the Sudakov exponent is quite 
simple and at the same time it incorporates soft wide-angle gluon emissions exactly. 

Consider next the reweighting factor W, whose main purpose is to provide a Sudakov 
suppression factor multiplying the matrix element according to a synthetic shower that 
corresponds to the Born level final state {p, f} m . This factor is not present in either of the 
existing programs il and [|, |TcJ that add parton showers to an NLO calculation. The 
reweighting factor W is simply taken from the paper of Catani, Krauss, Kuhn and Webber 
|TT| 1 and adapted to the present case of an NLO calculation. There may be improvements 
possible. For instance, Lonnblad 0 and Mrenna and Richardson [jl8| have provided more 
sophisticated algorithms for calculating a function that plays this role but more closely 
matches the Sudakov factor that would actually be generated in a shower Monte Carlo 
program. 

Consider finally the factor / that represents the expectation value of the observable 
after secondary showering and hadronization that is initiated from the partons produced at 
the Born level together with one hardest splitting. We have left this function as something 
that can be freely chosen subject only to a few very general conditions. If one wants a 
HERWIG style angle-ordered shower, the vetoed shower treatment of Catani, Krauss, Kuhn 
and Webber [ 111 is useful. Additionally, the recent paper of Nason [15] explains in some 
detail how an angle ordered-shower can be initiated after a first step consisting of the 
hardest splitting rather than the widest angle splitting. If one would like a Pt ordered 
shower that properly incorporates coherence, the recent paper of Sjostrand and Skands 
114] describes such a treatment, including the feature that emissions are generated in the 
order of decreasing overall Pt in the whole event rather than Pt along a single line. If one 
would like a shower based on the Color-Dipole Cascade Model of ARIADNE, then the work 
of Lonnblad [17] is to be recommended. We note that the most natural choice of shower 
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to follow the initial splitting based on the Catani-Seymour dipole splitting functions and 
kinematic definitions would be a shower in which all of the splittings are based on this 
choice. A construction based on this choice is, however, beyond the scope of this paper. 
Instead, we leave the choice of secondary showering scheme open. We emphasize that the 
same NLO calculation with the hardest splitting included could then be coupled to shower 
Monte Carlo programs using a variety of splitting styles as long as the subsequent shower 
is initiated so as to account for the splitting that has already occurred. 

The algorithms covered here include final state showers, but do not allow for strongly 
interacting partons in the initial state. The dipole subtraction scheme of Ref. |li| covers 
initial state partons as well, but we thought it best to leave the case of initial state showers 
for a separate publication. 

We certainly expect that the algorithm presented here can be improved. For instance, 
we note that the Air-jet matching scheme that we have used is rather like the phase space 
slicing method in NLO calculations and we wonder if it might be replaced by something 
that is more like the subtraction method. Additionally, it would be good to incorporate 
knowledge of one loop splitting functions into the algorithm so as to better match the one 
loop perturbative calculations. 

We have not constructed computer code that could demonstrate that the algorithm 
presented here works in practice. There may well be undesirable features that show up 
when we have real code to test. However, the algorithms presented here are a refined 
version of the algorithm presented in Refs. 0, pf . This cruder algorithm does exist as code 


and, when coupled appropriately to PYTHIA, appears to work quite well [1C]. 
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